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Abstract 

There is evidence that ice age cycles are paced by astro- 
nomical forcing, suggesting some kind of synchroniza- 
tion phenomenon. Here, we identify the type of such 
synchronization and explore systematically its unique- 
ness and robustness using a simple paleoclimate model 
akin to the Van der Pol relaxation oscillator and dy- 
namical system theory. As the insolation is quite a 
complex quasiperiodic signal, the traditional concepts 
of phase- or frequency-locking used to define synchro- 
nization to periodic forcing are inadequate. Instead, we 
explore a different concept of generalized synchroniza- 
tion in terms of (coexisting) synchronized solutions for 
the forced system, their basins of attraction and insta- 
bilities. 

We propose a clustering technique to compute the 
number of synchronized solutions, each of which cor- 
responds to a different paleoclimate history. In this 
way, we uncover multistable synchronization (reminis- 
cent of phase- or frequency-locking to individual peri- 
odic components of astronomical forcing) at low forc- 
ing strength, and monostable or unique synchronization 
at stronger forcing. In the multistable regime, different 
initial conditions may lead to different paleoclimate his- 
tories. To study their robustness, we analyze Lyapunov 
exponents that quantify the rate of convergence towards 
each synchronized solution (local stability), and basins 
of attraction that indicate critical levels of external per- 
turbations (global stability). We find that even though 
synchronized solutions are stable on a long term, there 



exist short episodes of desynchronization where nearby 
climate trajectories diverge temporarily (for about 50 
kyr). We also show that when a synchronized solution 
approaches the boundary of its basin of attraction, small 
perturbations may cause a jump to a different (coexist- 
ing) paeleoclimate history. 

Our study brings new insight into paleoclimate dy- 
namics and reveals a possibility for the climate sys- 
tem to wander throughout different climatic histories re- 
lated to preferential synchronization regimes on obliq- 
uity, precession or combinations of both, as environ- 
mental parameters varied throughout the history of the 
Pleistocene. 

Keywords : Climate Models Milankovitch Oscil- 
lator Generalized Synchronization Lyapunov expo- 
nent 



1 Introduction 

This article is a contribution to the field of pa- 
leoclimate dynamics theory, which has experienced 
many developments in terms of ice age models 
since many years notably by |Le Treut and Ghil, 1983| 
ISaltzman and Maasch, 199Q| , and many others, and re- 
mains an active research field. Paleoclimate modeling 
is a complex problem, hence an uncomfortable situa- 
tion for a scientist. On one hand, the data are scarce 
and marred by uncertainties. On the other hand, there is 
not a single well established model, the problem is non 
autonomous, the forcing is aperiodic, and stochastic ef- 
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(a) The LR04 stack | Lisiecki and Raymo, 2005) of 57 benthic 5^^0 (b) CO2 composite record [ppmv] |Luethi etal., 2008| . High values 
[%o] records; the is a proxy for the global volume of ice. High of CO2 correspond to a warmer climate (interglacial state), 
values of 5^^0 correspond to a colder climate (glacial state). 



Figure 1 : The long-term climatic signals reveal a slow-fast dynamics (only the most recent 500 kyr of the data are 
displayed here). The areas in grey (width Tsiow) are wider than the areas in white (width Tfast)'- while glaciation 
is a slow process of ice build-up, deglaciation occurs much more rapidly (t^/^w > ^fast)- One also recognizes the 
last deglaciation which started some 20 kyr ago, up to the present time (^ = 0). 



fects are present. 

Here, we focus on the slow variations of climate 
over the few last million years, which include the 
phenomenon of ice ages [Hays et al., 1976|, that is, 
the repeated growth and decay of ice sheets in the 
Northern Hemisphere of a total mass as big as mod- 
ern Antarctica's. When examining long-term cli- 
matic signals like the 5.3 Myr-long stack produced 
in [ Lisiecki and Raymo, 2005 1 , or the 800 kyr-long 
EPICA Dome C Ice Core from |Luethi et al. 72 008], 



plotted respectively in Fig. |l(a)| and Fig. |l(b)| for the 
last 500 kyr, one immediately identifies three clearly 
visible features of the climatic time series: 

(i) oscillations', the signal oscillates between higher 
and lower values of ice volume corresponding to 
the glacial and interglacial states. 



(ii) asymmetry: in Fig. 1 1(a) [ typical transitions from 
a minimum to a maximum take much longer 
than transitions from a maximum to a mini- 
mum: deglaciations occur much more rapidly 
{'^fast ^10 kyr) than glaciations (t^/^w ^ 80 kyr), 
giving a distinctive saw-tooth structure in the 
glacial/interglacial (G/I) cycles, especially pro- 
nounced over the last 500 kyrs. 



(iii) 100-kyr dominant period: this has 

been identified by many authors since 

proecker and van Donk, 19701 . Note that 
the G/I cycles are not periodic. 

The asymmetry in the oscillations has been studied 
by many authors. In order to reproduce it, some authors 
use underlying physical principles to build phenomeno- 
logical models that exhibit slow-fast dynamics reason- 
ably mimicking the climatic proxies |Saltzman, 2002| |. 
Others assume this asymmetry by explicitly defining 2 
different parameters such as the time intervals Tup = 
Tsiow and Tdown = ^fast | Ashkenazy, 2006| or time con- 
stants {Tr and Tf in tPaillard, 199 8 J and and Tc in 
llmbrie and Imbrie, 1980| ). Whatever the model, it has 
to ultimately exhibit asymmetric oscillations under the 
effect of the forcing, as it is aimed to mimic the oscil- 
lations between G/I states. Relaxation oscillators are 
therefore very straightforward natural candidates of ice 
age models. In this article, we will consider a slightly 
modified van der Pol oscillator model to illustrate the 
new contributions of our synchronization concepts. 

In this paper, we concentrate on the influence of the 
astronomical forcing on Earth's climate. This forcing 
is induced by the slow variations in the spatial and sea- 
sonal distributions of incoming solar radiation (insola- 
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tion) at the top of the atmosphere, associated with the 
slow variations of the Earth's astronomical elements: 
eccentricity (^), true solar longitude of the perihelion 
measured with respect to the moving vernal equinox 
(t(J), and Earth obliquity (e^). These quantities are now 
accurately known over several tens of millions of years 
iLaskar et al., 2QQ4| , but analytical approximations of 
e, ^sinCJ, and £e valid back to one million years have 
been known since | |Berger, 1978| . They take the form 
of d' Alembert series (£A/ sin[CL)/r + 0/]). 

The external forcing F{t) used throughout this article 
is the insolation at 65°N latitude on the day of the sum- 
mer solstice. That specific insolation quantity is com- 
monly related to the Milankovitch theory and can be 
thought of as a measure of how much ice may melt over 
summer. It can be written under the following compact 
form: 

35 

^(0 = Y^l^i sin((O,-0 + Q cos((0,-0] (1) 

i=\ 

where the value of the 3 x 35 parameters (including 
(Oi^Si^ and q) are given in the Table [T] of Appendix [A| 
The coefficients were extracted from | [Berger, 1978| by 
performing a linear regression of the insolation on the 
(Oi. The validity range of this approximation is [-1 Myr, 
Myr], and its mean error (mean of the absolute value 
of the difference, compared to [ Laskar et al., 2004 ]) is 
6.7 W/m^ with peaks at 27.5 W/m^. Note that the mean 
value (494.2447 W/m^) has been removed; the theoret- 
ical framework that allows to work with anomalies was 
justified by [ S altzman and Maasch, 1991 ]. In short, this 
can be done as we are interested in oscillations, and not 
in the mean values themselves. The quasiperiodi^nsi- 
ture of the insolation forcing F{t) is illustrated in its 
spectrum decomposition in Fig. |2] Precession is domi- 
nated by two harmonics around 19 and 23 kyrs (1 kyr = 
1 ,000 years) and obliquity is dominated by an harmonic 
with a period of 41 kyrs but it bears periods as long as 
1,200 kyrs. 

Synchronization 

There is ample evidence that the astronomical forc- 
ing influences the climate system. The phrase 'pace- 
maker of ice ages' was coined in a seminal paper 
[Ha ys et al., 1976| to express the idea that the timing 
of ice ages is controlled by the astronomical forcing, 

^ A quasiperiodic signal is the superposition of several periodic 
signals with uncommensurate periods. 



solstice at 65°N 



I Obliquity terms 
I Precession terms 



Figure 2: Spectrum decomposition of the incoming so- 
lar radiation (insolation) at 65°N latitude on the day of 
the summer solstice {F(t) given in Eq. This IS 

a graphical representation of Table [T] in Appendix [A| 
with aj = sj -\- cj and Ti = In/cOi. The eight largest 
parameters a/, which represent already 80% of the sig- 
nal, are the major components of the insolation; they 
clearly come from the precession (19 and 23 kyr), and 
from the obliquity (41 kyr) associated series. This inso- 
lation is the forcing used to construct all Figures subject 
to the quasiperiodic astronomical forcing. The main 
harmonic £\ associated with obliquity has an angular 
velocity (Og^ = 0.1532 rad/kyr (Je^ = 41.0 kyr) and an 
amplitude of a^^ = 11.77 W/m^. The three main har- 
monics associated with precession are denoted p\, pi 
and /73. 
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while the ice age cycle itself is shaped by internal sys- 
tem dynamics. The paradigm has prevailed since then 
and is it still supported by the most recent analyses 



need not be unique. If there are two or more rela- 
tionships ^ for the same parameter settings, we speak 



of palaeoclimate records [Lisiecki and Raymo, 2007 
|Huybers, 2007| . The notion of 'pacemaker' natu- 
rally evokes some sort of synchronization. How- 
ever, despite some attempts, the actual type of syn- 
chronization has not been clearly identified or demon- 
strated to date. For example, | |Ashkenazy, 2006[ 
Tziperman et al., 2006] speak of "nonlinear phase- 
locking" although they do not define suitable "phase 
variables" that can be used to demonstrate a fixed-in- 
time relationship between phases of the forcing and the 
oscillator response. 

Synchronization, as a universal nonlinear phe- 
nomenon, is a pervasive process in Nature, as it is asso- 
ciated with rhythmic processes. It is therefore not sur- 
prising to have synchronization also in Paleoclimatic 
Sciences. Depending on the forcing type (periodic, 
chaotic, stochastic), one can distinguish many types 
of synchronization including complete, lag, phase, 
frequency, identical, generalized [Rulko v et al., 19951 , 
achronal and isochronous [Wu et al., 20061 , and even 
noise synchronization. For a review, the reader is re- 
ferred to IBalanov et al., 20091 [Pikovsky et al., 2001 1 , 
among others. While some terminology is still debated, 
[Brown and Kocarev, 2000] proposed an unified defini- 
tion of synchronization for dynamical systems — there is 
synchronization if there exists a relationship h between 
the measured properties of the forcing, g{u), and those 
of the oscillator, giy): 



h{g{u),g{v))=0, 



(2) 



that is fixed-in-time, meaning that h is time inde- 
pendent. Because we are interested in synchroniza- 
tion that is stable, for arbitrary initial conditions u{0) 
and v(0) that do not satisfy we require that 
IBrown and Kocarev, 20001 : 



lim h{g{u{t)),g{v{t)))=0. 



(3) 



of multistable synchronization [Pikovsky et al., 2001 



Ch. 15.3.2]. Then, which of the relationships the system 
settles to will depend on initial conditions. 

In this paper, we use a simple van der Pol oscilla- 
tor model to identify and illustrate for the first time the 
phenomenon of generalized synchronization between 
ice age cycles and astronomical forcing. The dynam- 
ical systems approach outlined in the next section (i) al- 
lows for stability analysis of such synchronization, (ii) 
uncovers interesting effects relating to the robustness 
of the synchronization with respect to external pertur- 
bations, and (iii) uncovers the phenomenon of multi- 
stable synchronization that has been overlooked by pre- 
vious studies. We show that, in contrast to claims in 
[ [Tziperman et al., 2006| , synchronization needs not be 
unique. 

The article is structured as follows. Section [2] intro- 
duces a slightly modified version of the van der Pol os- 
cillator as a suitable model for studying synchroniza- 
tion of ice ages to astronomical forcing. In Section 
|3j we analyse synchronization to periodic forcing and 
quasiperiodic astronomical forcing in terms of largest 
Lyapunov exponents. Section|4]is dedicated to the study 
of multistable synchronization in terms of attracting tra- 
jectories in the phase space of the forced system, and 
the associated basins of attraction. In Section |5j we in- 
vestigate effects of the symmetry-breaking parameter j3 
for the van der Pol oscillator model. Section |6] is con- 
cerned with the robustness of the synchronization and 
focuses on two aspects relating to predictability. Firstly, 
it shows that the local stability can be lost temporar- 
ily causing divergence of nearby climatic trajectories. 
Secondly, it demonstrates that in the multistable regime 
external perturbations (such as noise) may cause jumps 
between coexisting synchronized solutions when these 
solutions come close to their basin boundary. To be 
clear, all the treatment below is deterministic, except 



for Figs. 14 and 16 



For example, if g{u) = u, g{v) = v, u and v have the 
same dimension, and ^ can be written as w = v, we 
speak of identical synchronization. More generally, if 
vectors u and v have different dimensions and ^ can- 
not be reduced to more than a functional relationship 
u = H{v), we speak of generalized synchronization; 
see also [ [Abarbanel et al., 19961 [Rulkov et al., 1 995', 
[Pikovsky et al., 2001| . Note that the relationship ^ 



This article requires some basics of Dynam- 
ical Systems theory (dynamical systems, nonlin- 
ear oscillations, limit cycles, bifurcations of vec- 
tor fields, etc.), for which we refer the reader 
to [ [Guckenheimer and Holmes, 1983[ [Arnold, 1983] 
[Strogatz, 19941 , to fSaltzman, 2002] for dynamical pa- 
leoclimatology, and also to [ [Savi, 20051 for a review of 
some useful concepts. 
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2 Generic ice age model: a modi- 
fled van der Pol relaxation oscil- 
lator 

The hypothesis at the basis of the work by Milankovitch 
IMilankovitch, 19411 is that changes in total amount of 
continental ice (say: x) are driven by summer insolation 
F{t) already described in Eq. ([T]). One straightforward 
interpretation of this hypothesis is a simple differential 
equation x = —cm^{x) /dx — yF{t), where d^^{x)/dx is 
the derivative of a climatic potential and yis the forcing 
efficiency. However, models of this form fail in practice 
to correctly capture the rapid deglaciation phenomenon. 
We therefore propose to model the paleoclimatic dy- 
namical system with a dissipative self-sustained oscil- 
lator resembling the classical van der Pol oscillatoij^ 

Ti = -[y^l5-YF{t)] (4a) 
ry = -a[^\y)-x] (4b) 

where ^'{y) =y^ /3—y. Note that this system is nonau- 
tonomous because the right-hand side depends explic- 
itly on time. 

The physical interpretation of the model is as 
follows. Ice volume x integrates the external forcing 
F{t) over time but with a drift j + j3. Assuming 
a 1, J is the faster variable whose dynamics is 
controlled by a two- well potential ^{y). For exam- 
ple, there are arguments that the dynamics of the 
Atlantic ocean circulation may be approximated by 
. 4b|| Rahmstorfetal., 20051 
Further interpretation and 
variable can be found in 



an equation similar to Eq 
Dijkstra et al., 2003| |. 
discussion of the fast 
I Saltzman et al., 1984| Tziperman and Gildor, 2003 



[Paillard and P arrenin, 20041 [Tziperman et al., 2006 



[Crucifix, 2011]. The parameter T sets the slow time 
scale. The coupled system Eq. ^ has one stable 
equilibrium solution for > 1 and a stable periodic 
orbit for < 1. The ratio of time spent near the 
two stable branches of the slow manifold given by 
^\y) = X depends on j3 (see the paragraph "Time 
Spent" in Appendix |B]). We use Tulc to denote the 
period of the stable periodic orbit and COulc = '^^/Tulc 
to denote the corresponding angular velocity. 

Relaxation oscillators have been pro- 
posed previously to study ice ages 
ISaltzman et al., 1984| [Tziperman and Gildor, 2QQ3| 



[Paillard and Parrenin, 2004| although, to our knowl- 
edge, in a less general form than here. We adopted 
this forn£] because it is very close to the well-studied 
van der Pol oscillator, and a good agreement (timing of 
glaciations and deglaciations, and their amplitude) with 
ice volume proxies was easily found for well chosen 
values of a,j3,7 and t (Fig. jsj. We note, though, that 
small changes in parameters or additive fluctuations 
may easily shift the timings of ice-age terminations for 
reasons that will be clarified later in the paper. 

The definition of synchronization can be applied to 
our model Eq. ^ as follows. The astronomical forc- 
ing F{t) corresponds to u{t), and the state vector whose 
two components are the slowly- varying ice volume x 
and the faster variable y corresponds tov{t). For nonpe- 
riodic forcing, relationship ^ can be very complicated 
(non-functional or even fractal-like) and hence difficult 
to detect. Therefore, other methods of detecting ^ had 
to be developed. As suggested by the auxiliary system 
approach [ Abarbanel et al., 1996| |, relationships ^ and 
^ are implied by an (invariant) attracting trajectory in 
the {x^y^t) phase space of the nonautonomous forced 
system ^ [ [Wieczorek, 201 1| . In the remainder of the 
paper, such an attracting trajectory is denoted with AT 
and referred to as an attracting climatic trajectory or 
synchronized solution. All other solutions to Eq. ^ 
will be referred to as climatic trajectories. 

Previous approaches to nonlinear dynamics of 
quasiperiodically forced oscillators focused on 
discrete-time mappings and two-frequency forcing 
iGlendinning and Wiersig, 1999 , Osinga et al., 2000 
[Belogortsev, 1992, ,Broer and Simo, 1 998 |. They 
uncovered interesting dynamics including Arnol'd 
or mode-locked tongues consisting of 'interlocking' 
bubbles and open regions of multistability, nonsmooth 
bifurcations, and strange nonchaotic attractors. Here, 
we consider quasiperiodic forcing with 35 frequency 
components and focus on the regions of mode locking. 
Our approach is based on instabilities of attracting 
trajectories in the (x^y^t) phase space of the continous- 
time forced system because they relate directly to the 
concept of generalized synchronization. We can pro- 
vide a systematic study of generalized synchronization 
to astronomical forcing by demonstrating existence of 
such trajectories and exploring their local and global 
stability properties. More specifically, we perform three 
types of calculations. Firstly, a clustering detection 



^ We give in the Appendix|B]a summary description of the classi- 
cal van der Pol oscillator and of its dynamical behaviour. 



^Note that the van der Pol oscillator model is also used as a ref- 
erence in ISaltzman, 2002| page 101] for a coupled ocean/sea-ice 
model. 
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Figure 3: Top: the insolation forcing F(t), scaled by 
in order to work dimensionless. Bottom: the x and y 
climatic trajectories obtained using system Eq. ^ with 
a = 11.11, P = 0.25, 7= 0.75 and t = 35.09. With 
these parameters, co^i = I.ScOjjlc^ where (Ogi is the an- 
gular velocity associated with the dominant harmonic 
of obliquity and (Oulc is the angular velocity associated 
with the unforced system's periodic orbit. Blue dots 
correspond to the Lisiecki and Raymo stack (LR04) de- 
scribed in Fig. |l(a)[ Time ^ = corresponds conven- 
tionally to the year 1950. The model is in good agree- 
ment with the ice volume proxy. 



technique uncovers parameter regions with monostable 
(unique) and multistable (non-unique) synchroniza- 
tion. Secondly, the largest Lyapunov exponent along 
AT quantifies its long- and short-term local (linear) 
stability. Thirdly, a basin of attraction of AT quantifies 
its global (nonlinear) stability. Finally, we remark that 
in the theory of nonautonomous dynamical systems, 
attracting trajectories in the {x^y^t) phase space are 
linked to a modern and more general concept of a 



pullback attractor [Kloeden, 2000, Langa et al., 2002 
IKosmidis and Pakdaman, 2003nWiggins, 2003| . 



3 Synchronization of the paleodi- 
matic system to the insolation 
forcing 

Illustration of the synchronization phenomenon 

A typical climatic trajectory for 7=0 is shown in 
Fig. |4j from two different points of view: the time se- 
ries and the phase space portrait. In the time series, we 
recognize the slow variable x (the ice volume), while 
y exhibits slow-fast dynamics. This climatic trajectory 
is also shown in the two-dimensional (x, y) phase space 
of the autonomous system where arrows indicate direc- 
tion of the flow. The trajectory converges to the limit 
cycle with slow-fast dynamics (the speed along the tra- 
jectory can be visually assessed by the circles of the 
residence plot). Let us now consider a set of 70 ran- 
dom initial conditions in the (x,j) -plane at time = 0, 
and study the resulting climatic trajectories in the three- 
dimensional phase space {x^y^t) of the nonautonomous 
system (Fig. |5(a)| ) for time t > t^. One clearly sees that 
all trajectories converge to a cylinder — the attracting set 
in the {x^y^t) space. 

However, if we consider now an external forc- 
ing (7 > 0) then synchronization onto this forcing 
may occur under certain conditions [Ashke nazy, 2006[ 
Tziperman et al., 2006) . According to our definition 



([2j|3]), synchronization is represented by an attracting 
climatic trajectory in the {x^y^t) phase space. 

Consider first the case of a purely periodic forcing 
with a period of 7> = 41 kyr and strength 7 = 3.33. 
The 70 initial conditions give rise to climatic trajec- 
tories that, after a sufficiently long integration time, 
converge to two attracting trajectories (see Fig. |5(c)| ). 
Both attracting trajectories are periodic with period of 
Tr = 27> = 82 kyr, and time-shifted versions of each 
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(a) Time series : x (the ice volume) is the slow variable, while _y ex- 
hibits slow-fast dynamics. The same color convention white/grey as 
in Fig. ^has been used in order to highlight the slow and fast episodes. 
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(b) Phase space portrait in the two-dimensional (x,_y) phase space of 
the autonomous system: dynamical flow, limit cycle and residence 
plot (circles are spaced at every 0.5 kyr). More time is spent along 
the slow manifold (red dashed curve, corresponding to the function 
<I>'(_y) = _y^/3 — _y = x). The trajectory converges to the limit cycle 
(thin curve) with slow-fast dynamics. 



Figure 4: Dynamics of the unforced ice ages model Eq. ^ with a = 11.11, j3= 0.25, 7=0 and T = 35.09. 
The slow-fast variable is y, while x is always slow. Example of a typical climatic trajectory with initial conditions 
(^-500,J-50o) = (-0.24, -0.27). 



other. This phenomenon is described in the litera- 
ture as 2 : 1 phase-locking ox frequency -locking. Gen- 
erally speaking sl n : m frequency-locking is defined 
as a fixed-in-time relation between the frequencies of 
the forcing (cOf) and the oscillator response (cOr) of 
the form n(OR = mcO/r, where m and n are integers 
IPikovskyetal., 2001| p.52]. 



Then consider the case of the quasiperiodic insola- 
tion forcing described in Eq. ([T]) with 7 = 0.75 and 
T = 43.86. Figure 5(e) shows that the 70 climatic trajec- 
tories now converge onto three attracting trajectories, 
which reveals that synchronization can be multistable 
IPikovskyetal., 2001| p.348], |Balanov et al., 2009| 
p. 94]. This phenomenon is described in the literature 
as mode-locking [ Svensson and Coombes, 200"9l . Note 
that because of the quasiperiodicity of the insolation 
forcing, these attracting trajectories are no longer peri- 
odic nor time- shifted versions of each other. The num- 
ber of attracting trajectories depends on many factors 
including the dynamics of the unforced system, the na- 
ture of the forcing F{t), and the amplitude 7 of the forc- 
ing. We will study this in more details in Sec.|4] 



Detection of synchronization by the way of the 
largest Lyapunov exponent (LLE or Xmax) 

Local or linear stability of an attracting cli- 
matic trajectory can be quantified with the largest 
Lyapunov exponent (LLE) denoted here as Xmax 
IBenettin et al., 1980| . The quantity Xmax is a measure 
of the (average) exponential rate of divergence {Xmax > 
0) or convergence (A^ax < 0) of nearby climatic trajec- 
tories. Therefore, a negative value of ^rnax indicates a 
locally attracting climatic trajectory or generalized syn- 
chronization IPikovsky et al., 2001] [Wieczorek, 2009| . 
A transition from X^ax < to X^ax = indicates a bi- 
furcation where the attracting climatic trajectory disap- 
pears and generalized synchronization is lost. Null and 
positive values of Xmax indicate lack of synchrony (pos- 
itive Xmax indicates chaos but this regime is not encoun- 
tered here). In the case of periodic forcing, computa- 
tions of Xfnax can be easily validated with more precise 
and reliable numerical bifurcation continuation tech- 
niques (see § '41 kyr periodic forcing' below). 



Long-term Xmax and short-term A^^^ LLE's 
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(e) with the quasiperiodic insolation (7= 0.75, T = 43.86 Tulc ~ (f) Section of Fig. |5(e)| at t = 550 kyr: three clusters of trajectories 

125 kyr) given by Eq. {T}: the trajectories converge now to three at- are identified, corresponding to the three attracting trajectories born 

tracting trajectories for a long time, revealing a multistable synchro- from multistable synchronization. 



nization. 



Figure 5: Illustration of three different paleoclimate dynamical regimes: (top) without any forcing, (middle) 
frequency-locking 2:1 onto a 41 kyr periodic forcing (symbol '+' in Figs. 6(a) and 6(c)[), (bo ttom) generalized 
multistable svncbroniyation on the nnasineriorlic insolation forcing rsvmbol ' x ' in Fips The ice 
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The largest Lyapunov exponent 
cally definecf]as fcOtt, 2002i | : 



is mathemati- 



lim 

|5Z(0)K0 



lim -m^; — 
t^oo t |5Z(0)| 



(5) 



where SZ= [5x^5y] are vanishing perturbations about 
X and y, respectively, governed by the linearization of 
system Eq. Whereas this classical Xmax is de- 

fined in long term limit ^ oo), one can also define 
lAbarbanel et al., 19911 a short-term version, A^^^, by 
considering a finite time interval H (H = 50 kyr will be 
considered in this article): 



1,H ^ 



lim — m ' ^ ' , ' 
|5Z(0)K0 H |5Z(0)| 



(6) 



While gives the average or long-term stability 
information, ^^^^ can tell us about the behaviour of 
nearby trajectories within a short time interval H. For 
example, X^nax < does not necessarily imply < 
for some suitably chosen H. The definition ^ will be 
useful in studying the robustness of generalized syn- 
chronization in Sec.[6l 

For computing X^ax^ as the system Eq. ^ and the 
Jacobian have an analytical form, tangent space meth- 
ods IKantz and Schreiber, 2004| can be used; technical 
details are given in the Appendix [C] 

Influence of the parameters 7 and Tulc 

The two particular types of synchronization illus- 
trated in Fig. |5(c)| and |5(e)| have been obtained for a 
fixed value of the amplitude 7 of the external forcing 
and of the natural period of the unforced paleoclimatic 
system Tulc- Now, we are equipped to achieve a much 
broader view of the dynamics by performing a paramet- 
ric study on these two parameters. The quantity plotted 
in Figs. |6(a) and 6(b) is the largest Lyapunov Expo- 
nent at 3 Myr, far from the transient behaviour so that 



Xjnt 



H=3 Myr 



may already been considered as a good ap- 



proximation of 



41 kyr periodic forcing 

^Even if differential versions of the LLE have sometimes been 
developed mainly for computational efficiency purposes, we however 
preferred within this article to stick on the original definition of the 
LLE, because it is more standard and there is no insistent need for 
lowering computation time in the present framework, as the number 
of degrees of freedom of the system is reduced. 



Fig. 6(a) I corresponds to the case of the 41 kyr pe- 
riodic forcing (Tp =41 kyr). The synchronization re- 
gion i^max < 0) is composed of several V-shape regions, 
called ArnoVd tongues (phase- or frequency-locking), 
originating at 1, 2, 3, etc. times the forcing period 7>. 
These regions correspond to 1:1, 2:1, 3:1 frequency- 
locking zones (3:2 and 5:2 can also be guessed). Pe- 
riodic solutions are found within these regions which 
originate generally speaking at Tuix: = {m/n) 7>. No 
synchronization is possible when 7 is zero but synchro- 
nization may occur already for infinitesimally small 
7. Then, for increasing 7, the synchronization re- 
gion widens and synchronization becomes more sta- 
ble up to an optimum value 7* of the forcing. When 
7 > 7*, the synchronization becomes less and less ef- 
fective, because at large 7 the system is too much 
steered away from its natural dynamics; it may even 
be driven into chaos at yet higher forcing amplitude 
IMettin et al., 1993| but this case is beyond our focus. 

In order to perform an accurate validation of the 
synchronization region given by the LLE {Xmax < 
0) method, we computed the main Arnol'd tongues 
boundaries with the more accurate numerical contin- 
uation methods such as AUTO | Doedel et al., 20091 . 
The case of periodic forcing [F{t) = sm{cot)] with 
j8 = has been already extensively studied in the 
literature, analytically assuming some approximations 
IGuckenheimer and Holmes, 1983| p. 70-75], and using 
numerical algorithms for pseudo-arc length continua- 
tion IMettin et al., 1993 |. The usual approach extends 
the original nonautonomous system by additional dif- 
ferential equations for the forcing so that the system 
becomes autonomous, and then explores the (ft), 7) pa- 
rameter space. Note that the asymmetry introduced here 
with the parameter j8 adds slightly more complexity 
and induces additional features to the diagrams doc- 
umented in these papers. In this way, we computed 
Arnol'd tongue boundaries as saddle-node of limit cy- 
cle bifurcations for the extended system with a = 1 1 . 1 1 
and P = 0.25. 

Superposition of LLE calculations and bifurcation 



boundaries in Fig. |6(a)| shows that the synchronization 
regions obtained with two different techniques match 
perfectly. This is a confirmation that the method based 
on the LLE works fine and we will be able to use it for 
the case of the quasiperiodic insolation forcing. Note 



that bifurcation boundaries are also drawn in Fig. |6(c) 
in order to stress the correspondence with yet another 
method of detecting synchronization that will be dis- 
cussed in Sec.m 
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Astronomical quasiperiodic forcing 

For the case of the quasiperiodic insolation forcing 



(Fig. 6(b)), the region of synchronization appears to 
be in one single piece with some indications of well- 
separated tongues (mode locking) at small 7. In other 
words, whatever the value of the natural period Tjjlc 
of the paleoclimatic system, it has a higher probabil- 
ity of being synchronized onto the insolation forcing. 
Of course, for very low values of 7, there is still no 
synchronization. Note that 7 has been downscaled by 
^insoi ^ 5 (compared to a sine = a/2/2 for the 41 kyr pe- 
riodic forcing), in order to keep a realistic range com- 
parison. 



4 Non uniqueness: multistability 
and basins of attraction 

The detection of synchronization using the LLE 
{^max < 0) gives only an Yes/No-type of informa- 
tion, without making any distinction between different 
tongues as this would require information about multi- 



stability. For example. Fig. 6(b) indicates synchroniza- 
tion for the parameter settings marked with the symbol 
' x" but gives no information about the number of at- 
tracting trajectories (we know that there are three dif- 



ferent attracting trajectories from Fig. |5(e)| ). To explore 
the problem of multistable synchronization, we propose 
a clustering method that not only allows us to detect 
synchronization, but additionally provides information 
about the number of attracting trajectories denoted here 
with N. 

Multistability Analysis: numerical estimate of the 
number of attracting trajectories N by a clustering 
technique 

Consider the case of the quasiperiodic insolation 
forcing with the three attracting trajectories, i.e., = 3 



(Fig. |5(e)| ). Although N can often be easily assessed vi- 
sually, we want to automatically detect and count the 
number of ATs. As a matter of fact, can be eas- 
ily estimated in the following way. Fix a time t that 
defines a two-dimensional (x, j)-section in the (x^y^t) 
phase space. Then start with a grid of initial conditions 
at some time to < t and take t — to sufficiently large so 
that all the initial conditions converge to the attracting 
trajectories at t. Since each AT is represented by a point 
on the (x,j) -section, the problem of counting attract- 



ing trajectories reduces to a simple clustering problem. 
We designed a suitable automatic cluster detection al- 
gorithm that counts the number of clu sters to obtain an 
estimate of N. For example. Fig. 5(f) show^the {x^y)- 
section of the three-dimensional {x^y^t) phase space at 
^ = 550 kyr, given 70 initial conditions at = 0. The 70 
trajectories converge onto three (highly concentrated) 
clusters corresponding to the three attracting trajecto- 
ries. 

The idea of using clustering analysis for paleocli- 
matic dynamics comes from the natural fact that clus- 
tering is another way of looking at generalized syn- 
chronization where negative LLE makes the trajecto- 
ries cluster more efficiently. This provides another in- 
sightful viewpoint on the problem of identification of 
number of synchronized solutions of the paleoclimatic 
system: the more stable the synchronization, the more 
efficient formation of clusters. 

Two important aspects of cluster analysi^have to be 
considered to avoid risks of mis-identification of clus- 
ters: 

• the notion of a cluster is based on the threshold 
distance Jj^Jthat has to be carefully chosen. If dr 
is chosen too large, there will be just one cluster 
including all points; if it is too small, no clusters 
will form with more than one point. 

• in order to have sufficiently well formed clusters, 
the time interval t — to must be chosen large enough 
so that the transient behaviour is gone; an illustra- 
tion of the convergence is given in Fig.|7] 

Depending on the type and amplitude of the forc- 
ing 7, we can have potentially a whole range of pos- 
sible number of attracting trajectories N, ranging from 
one |Tziperman et al., 2006| , to a few (two in the 41 
kyr periodic forcing example in Figs. |5(c)] and |5(d)( or 



three in the quasiperiodic insolation forcing example in 
Figs. |5(e)| and |5(f)|). When no forcing is considered 



^See also Fig.jlsjfor a detailed view. 

^Cluster analysis or clustering is the assignment of a set of obser- 
vations into subsets (called clusters) so that observations in the same 
cluster are similar in some sense. This is a common technique for 
statistical data analysis used in many fields for countless applications. 
There exists many types of clustering, along with several methods, 
among which: hierarchical, partitional, spectral, kernel PCA (princi- 
pal component analysis), k-means, c-means and QT clustering algo- 
rithms. 

^This threshold distance dr appears in any computation related to 
clustering analysis, or analogically Recurrence Plots (RP) analysis in 
complex networks, for determining neighbours |Marwan et al., 2009] 
[Donges etal., 2009| . 
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(a) Largest Lyapunov Exponent Xmax [kyr~^] for the 41 kyr periodic 
forcing case. The region with Xmax < corresponds to synchroniza- 
tion; we recognize its underlying Arnol'd tongue structure. The bi- 
furcation boundaries of the Arnol'd tongues obtained with the more 
accurate numerical continuation method AUTO are superimposed, for 
validation purposes (white curve), and match perfectly. 




40 60 80 100 120 140 
TuLC [kyr] 



(c) Numerical estimate of the number of attracting trajectories N for 
the 41 kyr periodic forcing case. In practice, small positive values cor- 
respond to synchronization onto a few attracting trajectories, while 
high values indicate no synchronization. Now the region inside the 
synchronization tongues is colored in function of N. For the tongue 
corresponding to a frequency-locking n: 1, we have n attracting tra- 
jectories. The bifurcation boundaries of the Arnol'd tongues obtained 
with the more accurate numerical continuation method AUTO are su- 
perimposed, for validation purposes (black curve), and match per- 
fectly. 




(b) Largest Lyapunov Exponent Xmax [kyr~^] in the case of the 
quasiperiodic insolation forcing given by Eq. {T}: the broad region of 
synchronization is a typical signature of the quasiperiodic insolation. 
The region of synchronization appears to be in one single piece with 
some indications of well-separated tongues (mode locking) at small 

r- 
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(d) Numerical estimate of the number of attracting trajectories N in 
the case of the quasiperiodic insolation forcing given by Eq. the 
structure is now much more complex, consisting of intermingled se- 
ries of Arnol'd tongues. The region with one attracting trajectory, 
corresponding to unique or monostable generalized synchronization, 
is the largest. However, there are also parameter sets with A/^ = 2, 3 or 
even more attracting trajectories. 



Figure 6: Detection of synchronization by two different methods: the largest Lyapunov Exponent {Xmax < 0) (top) 
and the numerical estimate of the number of attracting trajectories N by 3. clustering technique (bottom), and for 
two different types of forcing: a 41 kyr purely periodic forcing (left) and the quasiperiodic insolation forcing 

show for which values or the 
the climate system occurs for the ice age model Eq. ^ with a = 11.11, j3= 0.25, and T = 35.09. The symbol 
(+) refers to the specific climatic attracting trajectories illustrated in Fig. |5(c)| in the 41 kyr periodic forcing case, 
for which N = 2. The symbol (x) refers to the specific climatic attracting trajectories illustrated in Fig. 5(e) in 
the quasiperiodic insolation forcing case, for which N = 3. 
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Figure 7: Illustration of the importance of choosing the 
time interval t — t{) large enough so that the transient 
behaviour is gone, in order to have sufficiently well 
formed clusters. Here, choosing r = 550 kyr ensures 
that the eight clusters are already formed, starting from 
^0 = 0. 



Forcing 
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(Figs. 5(a) and 5(b) ), or there is forcing but no synchro- 
nization occurs, we find no clusters at all. This means 
that there are as many points in the (x,};) -section at time 
t as initial conditions at time ^o- Clearly, it is difficult 
to numerically distinguish between no synchronization 
and a large number of attracting trajectories (N ^ 1). 
Therefore, we restrict ourselves to just six different re- 
gions in Fig. |6] where we use white to indicate when 
there are none or more than five attracting trajectories. 
Now, we apply the numerical clustering analysis 



in the case of the periodic forcing (Fig. 6(c)) and of 
the quasiperiodic forcing (Fig. |6(d)|). We set r = 



and consider a grid of 49 initial conditions covering 
X e [-2.2; 2.2] and y e [-2.2; 2.2] at the initial time 
^0 = —40Tf for the periodic forcing (Tp is the period 
of the forcing), and = — 1600 kyr for the astronomical 
forcing. Two points in the (x,j) -section are estimated to 
belong to a different cluster if their Euclidean distance 
is greater than 0. 1 . 



41 kyr periodic forcing 

An illustration of the three possible synchronized so- 
lutions (N = 3) existing for the 3:1 frequency-locking 
on a periodic forcing is given in Fig. [8j where the 
response can be locked on one of the periods of the 
forcing. More generally, N corresponds to the num- 



Figure 8: Illustration of the three possible synchronized 
solutions (N = 3) existing for the frequency-locking 3:1 
on a purely periodic forcing, in the time series format. 
The response can be locked on one of the periods of the 
forcing. 



ber of forcing cycles associated with the synchroniza- 
tion regime (A^ = 1 for 1:1; = 2 for 2:1; = 3 for 
3:2, 3:1, etc.|j The resulting pattern of different N is, 
as expected, in agreement with the bifurcation diagram 
(Fig. |6(c)l ). For example, A/^ = 3 in the 3:1 tongue. This 
method allows one to visualize the 4:3 (N = 4) and even 
the 5:4 (N = 5) tongues to the left of the 3:2 tongue. It is 
also seen that N is generally larger where different syn- 
chronization regimes co-exist; this is the case between 
the 2: 1 and 3:1 regimes around 7=4. 

Astronomical quasiperiodic forcing 



Fig. 6(d) shows that synchronization occurs for 
most parameter configurations. The region with 
one attracting trajectory (N = 1), corresponding to 



^This statement relies on the system invariance with respect to a 
time- shift of one forcing period ( [Tziperman et al., 2006| show a very 
nice illustration of this point). 
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Figure 9: A careful examination of the lower part (low 



7) of Fig. 6(d) allows an insightful physical interpreta- 
tion. The underlying structure of the intermingled series 
of Arnol'd tongues can be understood: it is a mixing of 
several Arnol'd tongues series corresponding to each of 
the main components of the spectrum of the astronomi- 
cal forcing, yielding a resulting pattern which looks like 
misaligned combs. Subscales based on the period of the 
four main components (ei, p\, p2 and p^) of the insola- 
tion highlighted in Fig. [2] allows a clear recognition of 
each individual Arnol'd tongues series. 



unique or monostable generalized synchronization 
[Rulkov et al., 199 5], is the largest. However, there 
are also parameter sets with A/^ = 2, 3 or even more at- 
tracting trajectories. They indicate multistable gener- 
alized synchronization where different possible stable 
relationships ^ between the forcing and the oscillator 
response coexist. 

A closer view in the lower values of 7 is given in 
Fig. [9j which allows an insightful physical interpreta- 
tion. 

Three tongues with A/^ = 1 , 2, 3 are rooted at Tulc / Te^ 
= 1,2 and 3, respectively, suggesting a synchronization 
on the main obliquity component of the astronomical 
forcing of the same nature as synchronization on a peri- 
odic forcing. A series of other synchronization tongues 
with N > I appear; they correspond to 2:1 (A^ = 2), 3:1 
(A^ = 3), 4:1 (A^ = 4) and even 5:1 (N = 5) synchro- 



nization on the three leading components of precession, 
denoted pi, p2 and p^,. Consequently, the richness of 
the astronomical forcing effectively widens the param- 
eter range for which synchronization occurs, compared 
to a periodic forcing. The phenomenon may be under- 
stood intuitively: just as you are more likely to tune on 
some radio station if you are surrounded by a dozen of 
free FM emitters, the system is more likely to synchro- 
nize on the rich astronomical forcing than on a periodic 
forcing. Synchronization with A/^ = 1 or 2, found for 
larger 7, can be interpreted as a form of combined syn- 
chronization on both obliquity and precession. 

It is crucial to appreciate that synchronized solu- 
tions are not periodic and that, unlike the periodic forc- 
ing case, different synchronized solutions for a given 
set of parameters are not time- shifted versions of each 
other. The idea that different synchronized solutions co- 
exist is of practical relevance for paleoclimate theory. 
Namely, the set of parameters used to obtain the fit to 
the paleoclimatic records shown in Fig. [3] give two dis- 
tinct solutions at r = when started from a grid of initial 
conditions at = —700 kyr. Sensitivity studies show 
that the choice of ^ — is sometimes important for esti- 
mating correctly N. However, tests with t — to as large 
as 200 Myr of astronomical time suggest that several at- 
tracting trajectories may co-exist at the asymptotic limit 
of ^0 — 

A similar numerical clustering analysis plot for j3 = 
0.6 is shown in Fig. [TO] in order to give an idea of 
the effect of this parameter (a more detailed analysis 
is performed in Sect. [5]). The main conclusion about 
the multistability remains, but the particular values of 
N change, as the intermingled tongue series are differ- 
ent. 

Evolving shape of the basins of attraction 

Each AT]- (/ = 1 . . . A^) has its own basin of attractioi^ 
[Barnes and Grimshaw, 1997) 1, that is defined as the set 
of all initial conditions in the (x^y^t) phase space that 
converge to that ATf as time tends to infinity. For our 
nonautonomous system Eq. (|4]), we can study basins of 
attraction in the (x, 3; ) -section for different but fixed val- 
ues of initial time to, and observe how they vary with ^o- 
A given initial condition at time lies in the basin of 
attraction of ATf if it approaches ATf as time tends to 
infinity. Technical details about the computation of the 



^ A more formal definition of the basin of attraction for 
nonautonomous dynamical systems is given in I Kloeden, 2000| 
|Langa et aC2Q02j . 
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chosen initial condition. In the case of the periodic forc- 
ing, the two Ars are roughly equally likely for all as 
could be guessed from Fig. |5(c)[ However, this is not 



140 



: 0.6 in- 



Figure 10: Same as Fig. 6(d) but now with j3 
stead of j3 = 0.25. The particular values of change, as 
the pattern of intermingled tongue series is different, but 
the main conclusion about the multistability remains. 



basins of attraction by use of the specific classification 
algorithm developed (see Fig. [T8| are given in the Ap- 
pendix |D] Basins of attraction are of major importance 
because they provide the information about global or 
nonlinear stability of synchronization. If we care about 
predictability, basin boundaries indicate when a change 
in the attracting climatic history is likely. 

The evolving shape of the basins of attraction is 
shown in Figs. 



11 and [12) for the case of 41 kyr pe- 



riodic forcing and the quasiperiodic astronomical forc- 
ing, respectively. The evolution is shown as a comic 
strip, where each subfigure has an x— axis ranging from 
-1.5 to 1.5, and a axis ranging from -2.5 to 2.5 (the 
axis labels have been removed for a better readability). 

In the case of a periodic forcing (two basins), the pat- 
tern repeats itself periodically (compare the t^ — Q kyr 
to ^0 = 40 kyr, and to to = 80 kyr subfigures) in Fig.[TT] 
However, in the case of the quasiperiodic forcing (three 
basins), the pattern is much more intricate and seems 
not to repeat itself for the time horizon considered here. 

The ratio between the area of a basin of attraction 
and the considered area of the phase space can be inter- 
preted as a probability to converge to the correspond- 
ing attracting trajectory when starting from a randomly 



the case for the quasiperiodic forcing where the prob- 
ability to reach the same attracting trajectory may vary 
significantly in time. For example, the yellow basin is 
rather small at = kyr but becomes much larger at a 
later time = 90 kyr. 

In the multistable regime, if an ATf happens to lie suf- 
ficiently close to its basin boundary, then small pertur- 
bations could make the climate jump to another (coex- 
isting) AT^-^^, reducing predictability. This phenomenon 
is illustrated in Sec.[6l 



5 Influence of the symmetry- 
breaking parameter (5 

As the parameter j3 controls the asymmetry of the 
glaciation/deglaciation saw-tooth structure (a higher 
value of j3 leads to an enhanced asymmetry), it is useful 
to investigate its effect. We have already indicated in 



Fig. 10 that multistability depends on j3 in the case of 
the quasiperiodic insolation forcing. A more systematic 
approach encompassing the whole range of j3 is shown 
in Figs. |13(a)| and |13(b)[ for the case of 41 kyr peri- 
odic forcing and the quasiperiodic astronomical forcing 
case, respectively. 

First consider the 41 kyr periodic forcing (Fig. 13(a)| ). 
To understand this Figure, recall that the unforced oscil- 
lator (i.e., 7= 0) has a stable fixed point for |j3 1 > 1 and 
a stable limit cycle for |j3 1 < 1. 

The system responds almost linearly to the forcing 
when |j8| is sufficiently large. This explains regions of 
unique synchronization (N = I) where only one climate 
response is possible. The system becomes excitable 
when is just slightly greater than one. If the forc- 
ing is large enough it will excite oscillations. In this 
case, N is equal to the number of initial conditions if 
synchronization is lost, or to a smaller number if syn- 
chronization occurs. 

Consider now the interval — 1 < j3 < 1. For this, 
keep in mind (i) that the period of the unforced oscil- 
lation varies by almost a factor of two within the range 
< |j3 1 < 1, and (ii) that synchronization requires some 
relation between the period of the unforced oscillations 
and the forcing period. Consequently, synchronization 
on the periodic forcing occurs only for fairly narrow 
ranges of j3 that are symmetric around zero. The fig- 
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Figure 11: Evolving shape of the two basins of attraction for system Eq. with a = 11.11, j3 = 0.25, 7= 3.33 
and T = 35.09. Case of a 41 kyr periodic forcing : the pattern is 82 kyr periodic (frequency-locking 2:1). The 
function (y) = j^/3 — j = x, corresponding to the slow manifold, is also shown (dashed curve). The attracting 
trajectories (symbol) sometimes lie close to the boundary of their own basin of attraction. 




Figure 12: Evolving shape of the three basins of attraction for system Eq. (|4]), with a = ll.ll,j3= 0.25, 7 = 0. 12 
and T = 35.09. Case of the quasiperiodic insolation given by Eq. ([T]): the pattern is quasiperiodic, constituting the 
specific signature of the insolation. The function (y) = j^/3 — j = x, corresponding to the slow manifold, is 
also shown (dashed curve). The attracting trajectories (symbol) sometimes lie close to the boundary of their own 
basin of attraction. 
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ure reminds us of Arnol'd tongues. The main synchro- 
nization regimes detected here correspond to 4:1, 3:1 
and 5:2 frequency-locking. Outside these synchroniza- 
tion regimes, the system fails to converge to a suffi- 
ciently small set of attracting trajectories, meaning that 
the forcing is not as an efficient pacemaker. 

Finally, compare this situation with that obtained 
with the astronomical forcing (Fig. 13(b)| ). Synchro- 
nization now occurs in a larger area of the parameter 
space. Whereas the structure of the periodic forcing 
is preserved as long as the forcing amplitude is low 
enough, there is a much more richer and complex pat- 
tern of different N for larger y. This pattern emerges 
from the interaction with different harmonics and their 
beatnotes. 



6 Robustness of synchronization 

Robustness or reliability of synchronization can be 
studied in terms of two properties of an attracting cli- 
matic trajectory. Local stability analysis based on the 
short-term LLE (A^^^) provides information about the 
short-term local convergence towards the AT. For ex- 
ample, a temporary loss of local stability indicated by 
^max ^ ^ will cause a temporary loss of synchrony and 
divergence from the AT even though the trajectory is 
stable on average {Xmax < 0)- Global stability analysis 
based on the geometry of basins of attraction for differ- 
ent AJs provides information about the system response 
to external perturbations such as random fluctuations. 
For example, an external perturbation may push a cli- 
matic trajectory outside of the basin of attraction of the 
AT. Robustness and uniqueness of synchronization are 
closely linked in the sense that the global stability is a 
factor only when there are coexisting attracting trajecto- 
ries. Robustness is compromised most when a tempo- 
rary loss of local stability coalesces with a weakening 
of the global stability. We will now briefly discuss these 
two effects that could restrict the prediction horizon for 
the evolution of climatic trajectories. 

Temporary desynchronization via loss of local sta- 
bility 

Some additional experiments made in our paleocli- 
matic framework reveal another strange behaviour in 
the system Eq. (|4), that can be deduced from a careful 
inspection of Fig. 14] In the presence of small additive 
noise, we notice that nearby trajectories could diverge 




Figure 14: The temporary divergence of nearby cli- 
matic trajectories reveals short-term local instabilities 
(e.g. around ^ ^ 157 kyr). This illustration was obtained 
by considering a set of 50 random initial conditions at 
^0 = within xxy e [—1,1] x [—1, 1], and by adding 
some noise of small amplitude in the system at every 
time step in order to trigger the instabilities. The model 
used is Eq. ^ with a = 11.11, j8 = 0.25, 7= 0.033, 
and T = 33.33. 



for some time, like those around t ^ 151 kyr. 

Such temporary divergence is similar to desynchro- 
nization bursts IRulkov et al., 19951 and strongly sug- 
gests to investigate the evolving sign of the short-term 
LLE A^^^ along the attracting climatic trajectory. We 
computed X^ax^^ along one of the two attracting tra- 
jectories of system Eq. (|4]), subject to insolation forcing 
given by Eq. ([T]). The result is shown in Fig. 15 where 
the attracting climatic trajectory has been coloured ac- 



cording to the values of 



H=50 kyr 



Although the sys- 



tem is synchronized on a long term (kmax = — 0.2 kyr~^ 
see Fig. [TT]), we see here that there exist episodes with 
positive values of the short- ternf^ LLE A^^^, revealing 
temporary desynchronizations [ Wi eczorek, 2009 | . This 
explains the divergence of nearby trajectories found in 
Fig.[T4l 

These results remained unchanged with respect to 
the most important parameters of the model. For ex- 
ample, our main conclusions about the stability re- 



For completeness, the very short-term {instantaneous) stabiHty 
has also been investigated (see Fig. [20|in the AppendixjEj, as a limit 
case H ^0, but it is less physically relevant within the paleoclimatic 
context. 
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Periodic forcing (41,000 years) astronomical forcing 




-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 

(a) (b) 



Figure 13: Numerical estimate of the number of attracting trajectories N plotted as a function of the asymmetry 
parameter j3 and the amplitude of the forcing y for the ice age model Eq. ^ with a = 11.11, j3 = 0.25, and 
T = 35.09, assuming (left) a 41 kyr periodic forcing and (right) the quasiperiodic astronomical forcing given by 
Eq. ([T]). Synchronization occurs in a larger area of the parameter space in response to the astronomical forcing 
than in response to the periodic forcing, and the pattern of different N for larger y is much more richer and 
complex; this structure emerges from the interaction with different harmonics of the astronomical spectrum and 
their beatnotes. 



Article submitted to Climate Dynamics 



B. De Saedeleer, M. Crucifix and S. Wieczorek 



Is the astronomical forcing a reliable and unique pacemaker for Climate? 



18 



-0.5 



I 



III 



main qualitatively valid, even for different values of 
a (like a = 100), or with a different type of potential 
(<J>^(j) = (j+1.7)(j+1.58)(3; + 0.8)(3;)(3;-0.5)),even 
if the shape and size of the limit cycle and the bound- 
aries of the basins of attraction are of course different. 
The effect of the insolation function F{t) has also been 
checked: we compared the attracting trajectories for the 
insolation given by Eq. ([T]), compared t o those for the 
insolation given by fLas kar et al., 2004| . As these in- 
solation functions are very similar, the results are also 
very similar, and no difference was noticed. 

At first glance, it may appear that these episodes of 
temporary divergence are not relevant to the robust- 
ness of synchronization because climatic trajectories 
converge back the the attracting trajectory on a long 
term. However, other effects may be present that could 
strongly amplify such temporary divergence. They are 
identified below. 

Sensitivity to perturbations: preliminary results 



I 



100 



120 



140 160 

t [kyr] 



180 



200 



Figure 15: Short-term largest Lyapunov Exponent 
^ax^^^^"^ of one of the two attracting trajectories of 
system Eq. (|4]), subject to insolation forcing given by 
Eq. ([T]), and for the same parameters as in Fig. 14 
The attracting climatic trajectory has been coloured ac- 
cording to the values of A^^^^^^^^, revealing tempo- 
rary desynchronizations; this explains the divergence of 
nearby trajectories found in Fig.[T4j e.g. around r ^ 157 
kyr. 



Consider again Fig. 12 showing -sections with 
coexisting attracting trajectories in the case of the 
quasiperiodic insolation, and their basins of attraction 
for different values of ^o- Suppose now that the system 
is subject to additive fluctuations (for example, these 
may represent volcanic eruptions). Under certain condi- 
tions, such external perturbations may cause a displace- 
ment of the trajectory to a different basin of attraction, 
causing a jumpp^to another attracting trajectory. 

As a further illustration of this idea we show in 
Fig. [16] two attracting trajectories (in the time series 



format) that coexist for the same system parameters as 
those used for the fit of Fig. [3] but with additive fluctu- 
ations added to the fast variable (see legend for details). 
A jump from one trajectory to another at around -475 
kyr (arrow) may clearly been identified. This shows 
that a climatic trajectory is robust against fluctuations 
if it stays away from the basin boundary but its robust- 
ness can weaken significantly due to the weakening of 
the global stability near the basin boundary. 

We conclude that externally triggered jumps between 
coexisting attracting climatic trajectories seem to be 
most likely when the temporary desynchronization due 
to the loss of local stability coalesces with the weaken- 
ing of the global stability due to the proximity to the 
basin boundary. 



In the periodic forcing case the phenomenon of jumping from 
one attracting trajectory to another in response to a perturbation is 
called a. phase slip [Pikovsky et al., 200T] p. 23 8]. 
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Figure 16: Sensitivity of the climatic system to per- 
turbations. The same solution of the ice ages model 
Eq. (|4]) as in Fig. [3] is plotted (black) along with one 
sample trajectory of the same system (red), but with 
additive fluctuations added to the fast variable: dy = 
-T"^(a(<J>'(j) -x))dt^b dW, with b = 0.5y^ and 
W a Wiener-process. The arrow shows the time of the 
jumping from one to the other climatic attracting trajec- 
tory, which reduces the predictability of the timing of 
the glaciations. 



7 Conclusions 

We have identified, illustrated, and provided a system- 
atic study of generalized and multistable synchroniza- 
tion between the climatic glacial/interglacial oscilla- 
tions and the astronomical forcing. For doing so, a se- 
ries of appropriate concepts and tools have been devel- 
oped. A van der Pol-type relaxation oscillator, designed 
to reproduce the slow-fast dynamics of the paleocli- 
matic records, has been used for illustration purposes, 
but the methodology proposed here may of course be 
applied to other paleoclimatic models. 

To study the uniqueness of synchronization, we pro- 
posed a convenient concept of the number of attracting 
trajectories in the phase space of the nonautonomous 
forced system, each of which corresponds to a syn- 
chronized solution. We computed the number of syn- 
chronized solutions using a numerical clustering tech- 
nique, and uncovered that in addition to a unique or 
monostable synchronization, there are parameter set- 
tings where one finds a nonunique or multistable syn- 
chronization. At low forcing amplitude we found re- 
gions of mode locking where the system synchronizes 
on the individual components of the astronomical forc- 
ing in a way that is similar to frequency-locking on pe- 
riodic forcing (Arnol'd tongues), giving rise to coex- 
isting synchronized solutions. As the forcing ampli- 
tude is increased, the combined effects of precession 
and obliquity restrict the number of possible synchro- 
nized solutions. The emerging stability diagram con- 
sists of a large region of monostable synchronization 
mixed with smaller regions of multistable synchroniza- 
tion. A comparison with periodic forcing shows that the 
system finds it easier to synchronize to quasiperiodic 
insolation forcing. It is therefore conceivable that the 
climate system wandered throughout preferential syn- 
chronization regimes on obliquity, precession, or com- 
binations of both, as environmental parameters varied 
throughout the history of the Pleistocene. 

The robustness of generalized synchronization was 
investigated in terms of the key indicators of stability of 
synchronized solutions: the long- and short-term largest 
Lyapunov exponent (local stability), and the geometry 
of the basins of attraction (global stability). We found 
that even though the synchronized solutions are locally 
stable on a long term, there exists episodes where the 
short-term largest Lyapunov exponent becomes posi- 
tive, leading to temporary desynchronizations. As a re- 
sult, climatic trajectories could diverge from the syn- 
chronized solution for some period of time (50 kyr typ- 
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ically). Moreover, we computed the evolving shape of 
the basins of attraction for the coexisting synchronized 
solutions, and uncovered that these solutions sometimes 
approach the basin boundary where they become very 
susceptible to external perturbations. As a result, a 
small perturbation could make the climate jump from 
one synchronized solution to another, reducing pre- 
dictability. Such jumps seem to be most likely when 
the temporary loss of the local stability coalesces with 
the proximity to the basin boundary. In this context, 
we briefly discussed the effect of stochastic perturba- 
tions on the timing of the deglaciations. We also illus- 
trated the difference between the evolving shape of the 
basins of attraction for periodic and quasiperiodic inso- 
lation forcing. In the case of the insolation forcing, we 
obtained an intricate pattern of basins of attraction that 
does not appear to repeat itself in time. 
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A Insolation model using 35 terms 

We give at the Table [T] the numerical values of the 3 x 
35 terms for computing the insolation F(t) following 
Eq. Qin Sect.[T] 



Table 1 : Coefficients for the insolation model described 
in Eq. ([T]) used for the incoming solar radiation anomaly 
at Summer Solstice and at 65°N. The first 15 terms cor- 
respond to the obliquity component, while the other 20 
terms correspond to the precession. 



(Oi [rad/kyr] 



Si [W/m^ 



Ci [W/m^ 



0.153249478547167 
0.158148666238883 
0.117190147169570 
0.155061775112933 
0.217333905941751 
0.150162587421217 
0.211709630908568 
0.156336369673117 
0.148350290855451 
0.206924898030688 
0.212525165090383 
0.229992875969202 
0.306498957094334 
0.311398144786051 
0.004899187691716 



-11.2287376815124 
-3.82499371467540 
2.28814805956066 
-1.29770081956440 
0.380973541305497 
1.54904176353302 
-0.810768209286259 
-0.918358442095885 
0.256895610735773 
-0.335783913402678 
0.267659228540196 
0.0696189733188958 
0.0247349748169616 
0.0138353727621181 
-0.160479848721994 



3.51682075211241 
-0.761851750263805 
1.80233702684623 
-0.635152963728496 
-1.46301711999210 

-0.0883941912769817 
-0.577980646565494 
0.196083726889428 
-0.524697312305024 

-0.0194792150128644 
0.128915417116900 
0.0746231714061285 
0.0140464395340974 
0.0304736668840422 
0.0594077968934257 



0.264933601588513 
0.280151350350945 
0.331110950251899 
0.328024059125949 
0.326211762560183 
0.269742342439881 
0.332923246817665 
0.371638925683567 
0.275366617473065 
0.323124871434233 
0.259396912994958 
0.324937167999999 
0.334197841377850 
0.274551083291250 
0.418183080135680 
0.111684123041346 
0.433400828898112 
0.126901871803777 
0.336010137943616 
0.177861471704732 



-15.5490493322904 
15.4319556361701 
9.0992249352734 
-7.87065384013669 
0.813786144754451 
0.0690448504314857 
1.44050770785967 
0.925324276580528 
0.997628846513796 
-0.378637986107629 
0.339477750517033 
-0.576082669762308 
0.346906064369828 
-0.441772417569753 
-0.0184884064645011 
-0.428006728186239 
-0.0049199219454561 
0.257509918217341 
-0.421809264016129 
-0.161827722328271 



-9.70406287110532 
4.75247271131525 
-10.6115244887390 
6.61544246063503 
-4.52641408099246 
-3.31639260969558 
1.06339286050120 
-1.02066758672154 

-0.362906496840039 
0.527217891742183 

-0.560509461538342 
1.18669572739338 

-0.648189701487285 
0.289576210423804 
0.109632390175297 
0.357006342316690 

-0.106148873639336 

-0.377639794223366 
0.324327509437558 

-0.362683869407858 



Appendix 
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B The classical van der Pol oscilla- 
tor : description and dynamical 
behaviour 

Classical van der Pol oscillator 

The classical van der Pol oscillator 
|van der Pol, 1926| is very well known, widely 
used and deeply studied; here the necessary description 
which is relevant to the scope of this article is given. 
For more details the reader is referred to the literature, 
see e.g. |Strogatz, 1994| . 

This van der Pol-type oscillator model can be re- 
garded as a special case of the FitzHugh-Nagumo 
(FHN) model |Kosmidis and Pakdaman7200 3|, also 
known as Bonhoeffer-van der Pol (BYP) model 
[Barnes and Grimshaw, 19971 . 

Historically, this model was derived by van der Pol, 
when he discovered the existence of oscillations in elec- 
trical circuits. He found that the oscillation period is de- 
termined by T* = RC (time constant of relaxation) in RC 
circuits, or by T* = in RL circuits ; hence he named 
this oscillation as relaxation oscillations. The interest- 
ing characteristics of the relaxation oscillation are the 
slow asymptotic behavior and the sudden discontinu- 
ous jump to another value. This oscillator has therefore 
been widely used in many fields like in physical and 
biological sciences, neurology, seismology, lasers, op- 
toelectronics, etc. 

This two-dimensional system exhibits the same basic 
dynamical features (limit cycle and slow-fast dynamics) 
necessary to fit the paleoclimatic data, except the asym- 
metry, for which the system will be corrected for, by the 
introduction of an additional parameter j8 (see Eq. 4a in 
Sect.[T]). 

This van der Pol oscillator is non conservative with a 
nonlinear damping, governed by the following second- 
order differential equation [Balanov et al., 2Q09| — in 
the case without forcing : 



:0 



(B.l) 



where /i > is a positive constant proportional to the 
damping. 

The van der Pol oscillator is a Lienard system 
IStrogatz, 1994| , because f{x) = — /i(l — x^) is an even 
function and because g{x) = x is an odd function (see 



Eq. B.l ), in the canonical form : 

x + /(x)i + g(x 



Moreover, since this Lienard system satisfies addi- 
tionally the Lienard theorem, it has a unique and stable 
limit cycl^^in the phase space surrounding the origin. 
This is a feature which is absolutely required in order to 
model the glacial-interglacial oscillations. 

By way of the Lienard transformation y = x — x^/?> — 



the second-order Eq. B.l can be transformed into 



an equivalent two-dimensional system of ODE's : 



y 



^{x-x^S-y) 
x/il 



(B.3) 
(B.4) 



which can also be expressed under the equivalent fol- 
lowing form (using the convention £ = — l//i, so that 
8^0 means more and more damping): 



X = {y-^\x))/£ 



introducing the function <l>'(x) : 



(B.5) 
(B.6) 



(B.7) 



the associated potential <J>(x) having the shape of a 2- 
well potential. By varying e, we can adjust the respec- 
tive time scales of the x and y variables; for e < 1, the 
slow-fast variable is x and the slow variable is y. Note 
that in the main body of this paper, x and y are inverted, 
so that X will represent the slowly varying ice volume 
(see Fig. [4]), as it is the classical convention for the dy- 
namical theory of paleoclimates. 

When /i > 0, the system exhibits a stable limit cycle, 
where energy is conserved. Near the origin x = x = 
0, the system is unstable and energy is gained, and far 
from the origin the system is damped and energy is lost 
(while when /i = 0, there is no damping, the solution 
is a pure harmonic signal, and there is conservation of 
energy everywhere). 

There is only one fixed point which is the origin 
(x,i) = (0,0). 

When X is small, f{x) ^ — /i (negative damping). 
Thus, the fixed point (x,i) = (0,0) is unstable (an un- 
stable focus when < /i < 2, and an unstable node oth- 
erwise); see IHilbom, 2000| . On the other hand, when 
X is large, f{x) ^ (positive damping). Therefore, 

the dynamics of the system is expected to be restricted 
in some area around the fixed point. 



(B.2) 



The nowadays still unsolved Hilbert's 16th Problem is related 
to determining the number and location of limit cycles for an au- 
tonomous planar vector field for which both functions are real poly- 
nomials of fixed degree. 
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When < /i <C 1 (small damping), the system can be 
rewritten in order to avoid division by jl. When jl^l 
(large damping), the oscillations become less and less 
symmetric. 

When driven, the van der Pol oscillator can lead to 
synchronization [Balanov et al., 2009] , but also to de- 
terministic chaos |Ruihong et al., 2Q(j8| , depending of 
the level of the driving force. Since the paleoclimatic 
system is driven by the insolation, one could also have 
synchronization and perhaps chaos - or at least some 
path on the routes leading to chaos; i.e. Lyapunov ex- 
ponents becoming positive. 

This oscillator, or slightly different versions 
of it (a similar one is the Poincare oscillator 
[Glass and Sun, 19941 ), has been mathemati- 
cally largely studied under many aspects: bi- 
furcation structure [Metti n et al., 1993| , fixed 
points and Arnol'd tongues, chaotic dynamics 
IChen and Chen, 2QQ8| |Parhtz and Lauterborn, 19871 , 
additive noise fDegli Esposti Boschi et al., 2002], 
basins of attraction [Ba rnes and Grimshaw, 1 997], 
etc. In this article, we focus on the synchronization, 
multistability and predictability properties of a slightly 
modified version of the classical van der Pol oscillator. 

But there is also a big difference of frameworks: usu- 
ally, a simple periodic forcing is considered, while in 
the case of the astronomical forcing, the forcing has a 
much more complex form (quasiperiodic), like the one 
given in Eq. ([T]). 

Time spent on the slow manifold 

In the case of large damping (/i 1, or £ ^ 0), the 
oscillations become less and less symmetric, and signif- 
icant differences becomes to appear between Tsiow and 
Tfast, as defined in Section [T] 

One can more precisely characterize the slow-fast dy- 



When the system enters the regioip] where \y 



^'{x) \ ^ we then have (|i| ~ e) < 



> £) : |i| an 



namics of the system by considering Eqs. B.5 - B.6 We 
see that we have in any case \y\ ^ £, which means that 
y is always slow. 

When the system is away from the curve y = <I>'(x), 
we have (|i| ~ l/e) {\y\ ~ e) : the vector field 
is mostly horizontal. And the travel speed is v = 
Y^'^M-^P ^ l/e ^ I, so that the system moves 
quickly in the horizontal direction (y is the slow vari- 
able, and X is the slow-fast variable). lfy> ^'(x), then 
i > 0, and the trajectory moves clockwise. See Fig. |4] 
for an illustration and [ Strogatz, 1994| for a discussion 
about the separation between the two time scales. 



\y \ are about the same order of magnitude, and the travel 
speed is V ~ e ^ : the system moves slowly along the 
curve y = <I>'(x), and eventually exits from this region, 
and so on : the system has a stable limit cycle. 
Hence we have the following relationships: 

Tfast^^/l^ or Tfast^e • (B.9) 

The period of the oscillations (or of the unperturbed 
limit cycle), given by 



Tcycle — '^slow ~^ '^fast 



(B.IO) 



is mainly determined by the time during which the sys- 
tem stays around the curve y = <I>'(x), which can be 
roughly estimated to be T oc 1/e, as v ~ e ^ 0. So, 
the larger the damping (lower e), the longer the period 
of the oscillations; that's why an additional time scaling 
(factor t) must sometimes be done in order to keep the 
100 kyr period for the limit cycle (see Fig [4]). 

Let us mention also that analytical expressions 
have been derived for the amplitude and period 
of the limit cycle |D'Acunto, 2006| in both cases 
({1 small or big), and also for the slow manifold 
equation |Ginoux and Rossetto, 20061 . There exists 
also a very fast transition upon variation of a pa- 
rameter j3, from a small amplitude limit cycle to 
a large amplitude relaxation cycle, explained by 
the so-called canard phenomenon, cycles, and ex- 
plosion IBenoit et al., 1981|[Guckenheimer et al., 20001 
[Guckenheimer and Haiduc, 20051 . 



C Technical details about the cal- 
culation of the Lyapunov expo- 
nents 

We first refer the reader to the seminal papers 
IShimada and Nagashima, 1979| [Benettin et al., 1980| . 
The methods for computing the Lyapunov (Char- 
acteristic) Exponents (LCE) vary depending on 
the fact that one wishes to achieve the full spec- 
trum of the LCE IWolf etal., 19851 , or only the 



^^This region is sometimes called "boundary layer" 
jGuckenheimer and Holmes, 1983| , by analogy with Fluid Dy- 
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largest one |Rosenstein et al., 1993| |. Analyti- 
cal derivations of the LCE's of the van der Pol 
oscillator do also exist |Grasman et al., 2005| . 
In this research, we computed Xmax using the 
standard method involving a Gram-Schmidt Re- 
orthonormalizaton (GSR) of the 'tangent vectors' 
[Shimada and Nag ashima, 1979| [Benettin et al., 19801 
[Wolf et al., 1985]; which is described in the review 
paper [ Ramasubr amanian and Sriram, 2000| . 

We remind here the fundamental principles. Con- 
sider an /t-dimensional continuous-time dynamical sys- 
tem : 

f=/(z,r) (CD 

where Z and / are ^-dimensional vector fields. To deter- 
mine the LCE's corresponding to some initial condition 
Z(0), we have to find the long term evolution of the axes 
of an infinitesimal sphere of states around Z(0). That is 
to say that we assume [ Qtt, 2002] 5Z{t) = 5Z(0) e^, 
with A = diag{X\ , . . . , A„), where Xi are the eigenvalues 
of the system. For this, we consider the linearization of 



Eq. C.l given by: 



d8Z 
dt 



J5Z 



(C.2) 



where / is the ^ x ^ Jacobian matrix defined by Jij = 
dfi/dZj. Then, starting from a uni t vector 5Z(0), the 
original system given by Eq. CI is integrated for Z 



together with the 5Z tangent system given by Eq. |C.2 
The evolution of 5Z is such that it tends to align with the 
most unstable direction (the most rapidly growing one). 
The choice of the initial vector of the tangent manifold 
may influence the convergence, but in practice a spin- 
up phase can be performed in order to find the good 
direction. 

The largest Lyapunov exponent ?^ax is then defined 
as : 

t \5zm 



lim lim 

|5Z(0)K0 t- 



It is of course impossible in practice to go to infin- 



it>l_} so the computation is always truncated to some 
finite final time, usually of the order of 10^~^ times the 
period of the forcing. The convergence can be more 
rapidly achieved if some transient behaviour is skipped 



(like illustrated in Fig. 17), i.e. we compute the LCE's 



^^The proof of the existence of such a Hmit has been made by 
lOseledec, 1968j . 



only when we are quite sure to be on the attracting tra- 
jectory. 

Note that there exists a whole spectrum of Lyapunov 
exponents A; (/ = 1,2, ...,^), that can be computed with 
a unit vector basis, and with a renormalization proce- 
dure. Although we are mostly interested in A^^^ in this 
article — a positive Amax is associated to a desynchro- 
nization — , our subroutine allows to compute all the 
spectrum of A^'s of any system. For the sake of flexi- 
bility, we used a symbolic software, so that the model 
could be very easily changed, and all functions (like the 
Jacobian matrix) are automatically derived once the ini- 
tial system is given. 

One of the standard and popular methods to com- 
pute the Lyapunov spectrum of a dynamical system in- 
volves a Gram-Schmidt Reorthonormalizaton (GSR) of 



the 'tangent vectors' [[Shimada and Nagashima, 1979 
IBenettin et al., 19 801 IWolf etal., 1985J; differential 
versions of this have also been formulated. Our subrou- 
tine includes the GSR procedure, which is required in 
order to avoid computational overflows, and degeneracy 
into a single vector. The frequency of reorthonormal- 
ization is not critical; as a rule of thumb, GSR is usu- 
ally performed on the order of once per characteristic 
period. Here, the normalization time step has been cho- 
sen optimally, that is to say the largest possible which 
preserves the accuracy. Since the GSR never affects the 
direction of the first vector in a system, this vector tends 
to seek out the direction in tangent space which is most 
rapidly growing. 

Our subroutine has been validated 
by comparing our LCE's to those of 
IRamasubramanian and Sriram, 200Q| for several 
systems (driven van der Pol, Lorenz '63, etc.); 
the order of the accuracy achieved is l%o for 
the Lorenz systenf^ One also checked that 
div / = Tr{^) = -a -1-13 = -13.66 = I^i A,- 
holds. 

Coming back to our system of Eqs 4a - 4b for which 
the Jacobian matrix is: 







1 



/z = ^{y,a) (C.3) 



we end up to the following results (Fig. [17]): the system 
Eq. ^ subject to the insolation (Eq. ([T])) is synchro- 
nized on a long term, since ^ax = —0.2 kyr~^ < 0. 



^^For the Lorenz system, one must pay attention to the 
spurious trivial set of LCE's corresponding to the origin 
[Shimada and Nagashima, 1979[|Bryant et al., r990j . 
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Figure 17: Convergence of the two LCE's of one of the 
two attracting trajectories of the system Eq. (|4]), subject 
to the insolation (Eq. ([T])); some transient is skipped. 
For T = 3.33, we have ?irnax = —0.2 kyr~^ < (for 
T = 35.09, the value must be scaled accordingly, which 
would give -0.019 kyr~^), so the system is synchronized 
on a long term. 



Some properties of the LCE's (AO 

The LCE's are very useful in order to characterize 
the dynamical behaviour of a system; for example, here 
are a few interesting properties, valid in the case of the 
dissipative system Eq. (pi in the autonomous case : 



• a dynamical system 
n has n LCE's and 



of 



dimension 
eigenvectors 



ILichtenberg and Lieberman, 1983| , 

• ^4=1 h = Tr{J) = div /, both being related to the 
growth of a volume of dimension n of the phase 
space IShimada and Nagashima, 1979| , 

• div(/) < for a dissipative system (hence the sys- 
tem has at least one negative exponent), 

• A; = along a limit cycle (tangent direction of the 
attractor) | Kantz and Schreiber, 20 04], 

• at least one LCE vanishes if the trajectory 
of an attractor does not contain a fixed point 
IHaken, 19831 , 

• the nature of the attractor can be described by an- 
alyzing the sign of the LCE's, which gives a qual- 
itative picture of the dynamics |Wolf et al., 1985 1 
[Miiller, 1995|| , also function of the dimension (ID 
to 4D), e.g. an attractor for a dissipative system 



with one or more A/ > is said to be strange or 
chaotic; if more than one A/ > 0, there is hyper- 
chaoticity |Rossler, 1979| . 



• the existence of A/ > is mathematically related 
to the theory of ergodicit}p^ of dynamical systems 
lEckmann and Ruelle, 1985| , 



• local bifurcations can be detected by detecting 
changes of signs of A^ . 

The LCE spectrum is closely related to the frac- 
tional dimension of the associated strange attractor by 
the Kaplan- Yorke formula [Kaplan and Yorke, 1979| |. 
There are a number of similar measures of the 
"strangeness" of strange attractors, like the frac- 
tal dimension |Theiler, 19901 , information dimen- 
sion, box-counting dimension, and the correla- 
tion dimension [Tsiganis et al., 1999^ and exponent 
|Grassberger and Procaccia, 1983], which allows to 
distinguish between deterministic chaos and random 
noise, and is computationally easier. 

It is also possible to estimate the LCE of a system 
by analyzing its time series, but with limited data, or a 
system subject to non negligible stochastic perturbation, 
the classical methods may provide incorrect or ambigu- 
ous results [ Ruelle, 1990], hence require specific meth- 
ods, like [ .McCaffrey et al., 1992||Liu et aL, 2005J . 



D Technical details about the com- 
putation of the basins of attrac- 
tion 

The practical computation of a basin of attraction is 
done as follows. 



Let us for example come back to the Fig. 5(e) with 
the three attracting trajectories ATi{x^y^t) due to the in- 
solation forcing. 

Now, we wonder which initial condition leads ulti- 
mately to which of the three ATf. This is the concept of 
the basin of attraction of a given A7^-, classically defined 
by the set of states that leads to a given ATf. Let us more 
precisely define the basin of attraction of a given AT] as 



^^An attractor is said ergodic (or transitive) if the points fill an en- 
tire physically realizable domain; intransitive if there are distinct at- 
tractors in several closed subdomains jSaltzman, 2002], and a system 
is almost intransitive if it resides for long periods of time in one or 
another domain that is not fully closed off, so that occasional exits 
from one domain to another can occur. 



Article submitted to Climate Dynamics 



B. De Saedeleer, M. Crucifix and S. Wieczorek 



Is the astronomical forcing a reliable and unique pacemaker for Climate? 



25 



the locus of all points in the (x, plane which lead to 
motion which ultimately converges on that AT]. 

We initialize many initial conditions on a fine rectan- 
gular grid covering the phase space. Each initial con- 
dition is then integrated forward to see which ATi its 
trajectory approached. If the trajectory approached a 
particular one of the three AT^'s, a dot coloured by the 
color identifying the AT] is plotted on the grid. For do- 
ing this, we need to define a target time at which we will 
do the classification, and a criteria for the classification. 

The classification algorithm is illustrated in Fig. [TSj 
where a cut has been made at a target time ^ = 550 kyr. 
The three AT] ' s are displayed, together with the location 
of the trajectories (black circles). To decide if a given 
trajectory ends up onto a given AT], we choose a max- 
imum distance from the AT] (dotted circle); taken here 
to be 1/4 of the minimum distance between two AT] 's. 
Trajectories falling into that dotted circle are classified 
as ending into AT] . 

We consider a given trajectory starting from = 
kyr, integrate it up to ^ = 550 kyr, and then we exam- 
ine its position with respect to the AT] at the same time 
T=550 kyr. If the distance to a given AJ^ is 'sufficiently 
small' (see the circles around the AT]), then this initial 
condition is coloured in the k^^ color, associated with 
the k^^ attracting trajectories AT]. 

Repeating this process for each initial condition on 
a fine grid of the whole phase space gives the shape of 
the basins of attraction. As we have three attracting tra- 
jectories, we will have three basins of attraction. Such 



basins of attraction obtained are shown in Fig. 19 



As we are in the case of a non-autonomous system, 
we then have to repeat this procedure for several starting 
times ^0 (the position of the AT] 's are constantly evolv- 
ing, hence the shape of the basins are also varying with 
time). This has been done to produce the evolving shape 



of the basins of attraction in Figs. 1 1 and 12 



Note that if t is too close from ^o, transient behaviours 
predominate, and not enough time has elapsed in order 
for the trajectory to be attracted by a given AT], hence 
the basins of attraction cannot be defined in that case. 

The glacial/interglacial cycles do exist since about 3 
million years, but only 8 limit cycles of 100 kyr period 
have been performed, since the Mid-Pleistocene Tran- 
sition (MPT). This should however be sufficient to con- 
verge onto the attracting trajectories, because the cli- 
matic trajectories are rapidly attracted on the limit cy- 
cle. So, if the ice age model Eq. (|4]) is a realistic one, it 
would be reasonable to state that the transient of climate 
dynamics has gone, and that we are currently probably 



• AT, 

AT, 

* AT, 



-1.5 -1 -0.5 



0.5 1 1.5 



Figure 18: Illustration of the classification algorithm, 
which allows to compute the basins of attractions. The 
position of the 70 trajectories (black circles) with re- 
spect to the three attracting trajectories (symbols) at 
time ^ = 550 kyr, starting from = 20 kyr, are shown. 
If the trajectory falls into a dotted circle, classification 
is possible. The 70 trajectories form three highly con- 
centrated clusters, corresponding to the three attracting 
trajectories. 




-2.5 -2 -1.5 -1 -0.5 0.5 1 1.5 2 2.5 

Figure 19: When integrating the system system Eq. (|4]), 
with many initial conditions (grid 201 x 121 = 24 321) 
starting from time = kyr, with the insolation forcing 
given by Eq. ([T]), we obtain the three basins of attraction 
above; each basin is coloured with the colour of the cor- 
responding AT] . The function ^' (y) = / —y = x, cor- 
responding to the slow manifold, is also shown (dashed 
curve). 
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somewhere on an attracting trajectory, if any. 



E Instantaneous stability 



As we considered the long-term (3 Myr) and short-term 
m = 50 kyr) local stabilities, we could wonder why not 
to consider also, at the other extreme, the 'very short- 
term' local stability (H ^ 0). Even if it has probably 
less relevance in the context of predictability of glacia- 
tions, it may however still be computed and understood 
as an instantaneous stability. As the horizon of time H 
tends to zero, it means that we no longer consider any 
motion in the phase space (no average in time), hence 
the instantaneous stability becomes also a local prop- 
erty in the phase space. 

The same formula (Eq. |6]) as for the short-term sta- 
bility is used, but now with an Horizon of // = 1 kyr, 
representing the very short-term. The result is plotted in 



Fig. 20 where the trajectory has been coloured with the 
A^^^^alues. There are black areas, which means local 
instantaneous instability. This black zone is always in 
the same region of the phase space: when y lies within 
the interval /* = [—1,1], whatever the initial condition. 
We now demonstrate mathematically why the black 



zone in Fig. 20 left is located within the interval /* = 

[-1,1]. 

Let us derive the theoretical expression of the exactly 
instantaneous LCE If we call A^ the eigenvalues 
of the jacobian matrix /, then the A/"'^^ are defined as 
^inst ^ li^^^Q i In A^ [ Drazinand Kiiigri992] , so that 
we can obtain the largest ^^^^^ by computing = 
max{3i{Ai)). 

As the Jacob ian m atrix / is a function of y only (and 
not x), see Eq. [C3 so do X^^^. The function XZ^{y) 
is plotted (Fig. |20[ right, black dotted curve), where 
we clearly see that it is positive on the interval [—1,1], 
which is precisely the same interval as /* . 





Figure 20: Very short-term stability. Left : the at- 
tracting climatic trajectory (initial conditions (xq, jo) = 
(—0.24, —0.27)) is coloured by the values of X^^^ 
the system Eq. Q with a = 11.11, j3 = 0.25, 7=0.033 
and T = 35.09. Instantaneous instability (^^^^ > 0) 
curs for }7 G /* = [-1, 1]. Right: the function 10 x X}^^^ 
(black dashed curve) has positive values for the same 
range j G /* = [— 1, 1]. For both graphs, the function 
^'{y) = X has been superimposed for an improved read- 
ing and easier interpretation. 



Moreover, the numerical values also match (e.g. the 
maximum of the function has a value of about 
0.28, which correspond to the maximum of the color 
scale). 

An interpretation of is that a particular trajec- 
tory makes a travel through the phase space 'painted' 
by X^ax^ integrating so the instantaneous stability to 
achieve the short-term and long-term stabilities. 



Article submitted to Climate Dynamics 



B. De Saedeleer, M. Crucifix and S. Wieczorek 



Is the astronomical forcing a reliable and unique pacemaker for Climate? 



27 



References 

[Abarbanel et al., 1991] Abarbanel, H. D. I., Brown, 
R., and Kennel, M. B. (1991). Variation of lyapunov 
exponents on a strange attractor. Journal of Nonlin- 
ear Science, 1(2):175-199. 

[Abarbanel et al., 1996] Abarbanel, H. D. I., Rulkov, 
N. R, and Sushchik, M. M. (1996). Generalized 
synchronization of chaos: The auxiliary system ap- 
proach. Phys. Rev. E, 53(5):4528^535. 

[Arnold, 1983] Arnold, V. (1983). Geometrical Meth- 
ods in the Theory of Ordinary Differential Equations. 
Springer- Verlag, New York, 1988 second edition. 
English translation of the original russian publica- 
tion: "Dopolnitel'nye Glavy Teorii Obyknovennykh 
Differentsial'nykh Uravnenii" (Additional Chapters 
to the Theory of Ordinary Differential Equations, 
Moscow: Nauka, 1978). 

[Ashkenazy, 2006] Ashkenazy, Y. (2006). The role of 
phase locking in a simple model for glacial dynam- 
ics. Climate Dynamics, 27(4):421^31. 

[Balanov et al., 2009] Balanov, A., Janson, N., Post- 
nov, D., and Sosnovtseva, O. (2009). Synchroniza- 
tion: From Simple to Complex. Springer- Verlag, 
Berlin, Germany. 

[Barnes and Grimshaw, 1997] Barnes, B. and 
Grimshaw, R. (1997). Analytical and numeri- 
cal studies of the bonhoeffer van der pol system. 
The ANZIAM Journal, 38(04):427^53. 

[Belogortsev, 1992] Belogortsev, A. B. (1992). 
Quasiperiodic resonance and bifurcations of tori in 
the weakly nonlinear duffing oscillator. Physica D: 
Nonlinear Phenomena, 59(4):417-429. 

[Benettinet al., 1980] Benettin, G., Galgani, L., 
GiorgilH, A., and Strelcyn, J.-M. (1980). Lyapunov 
characteristic exponents for smooth dynamical 
systems and for hamiltonian systems; a method 
for computing all of them, part 2: Numerical 
application. Meccanica, 15(l):21-30. 

[Benoit et al., 1981] Benoit, E., Callot, J., Diener, R, 
and Diener, M. (1981). Chasse au canard. Col- 
lectanea Mathematica, 31-32(l-3):37-l 19. 

[Berger, 1978] Berger, A. L. (1978). Long-term vari- 
ations of daily insolation and Quaternary climatic 
changes. /. Atmos. ScL, 35:2362-2367. 

Article submitted to Climate Dynamics 



[Broecker and van Donk, 1970] Broecker, W. S. and 
van Donk, J. (1970). Insolation changes, ice vol- 
umes, and the 018 record in deep-sea cores. Rev. 
Geophys.,Sil):169-m. 

[Broer and Simo, 1998] Broer, H. W. and Simo, C. 
(1998). Hill's equation with quasi-periodic forcing: 
resonance tongues, instability pockets and global 
phenomena. Soc. BrasilMat, (29):253293. 

[Brown and Kocarev, 2000] Brown, R. and Kocarev, L. 
(2000). A unifying definition of synchronization for 
dynamical systems. Chaos, 10(2): 344-349. 

[Bryant et al., 1990] Bryant, R, Brown, R., and Abar- 
banel, H. D. I. (1990). Lyapunov exponents from ob- 
served time series. Physical Review Letters, 65(13). 

[Chen and Chen, 2008] Chen, J.-H. and Chen, W.-C. 
(2008). Chaotic dynamics of the fractionally damped 
van der pol equation. Chaos, Solitons & Fractals, 
35(1):188-198. 

[Crucifix, 2011] Crucifix, M. (2011). Oscillators and 
relaxation phenomena in Pleistocene climate theory. 
Transactions of the Philosophical Transactions of 
the Royal Society ( Series A, Physical Mathematical 
and Engineering Sciences), In Press. 

[D'Acunto, 2006] D'Acunto, M. (2006). Determina- 
tion of limit cycles for a modified van der pol 
oscillator. Mechanics Research Communications, 
33(l):93-98. 

[De Saedeleer et al., 2010] De Saedeleer, B., Crucifix, 
M., and Wieczorek, S. (2010). Is the synchronization 
of the climatic system on the orbital forcing robust? 
Poster presented at the "SYNCLINE 2010: Synchro- 
nization in Complex Networks" conference, held on 
26-29 May 2010 at the Physikzentrum Bad Honnef 
(Germany). 

[Degh Esposti Boschi et al., 2002] DegH Es- 
posti Boschi, C, Ortega, G. J., and Louis, E. 
(2002). Discriminating dynamical from additive 
noise in the van der pol oscillator. Physica D: 
Nonlinear Phenomena, 1 7 1 ( 1 -2) : 8- 1 8 . 

[Dijkstra et al., 2003] Dijkstra, H. A., Weijer, W., and 
Neelin, J. D. (2003). Imperfections of the three- 
dimensional thermohaline circulation: hysteresis 
and unique-state regimes. J. Phys. Oceanogr, 
33:2796-2814. 

B. De Saedeleer, M. Crucifix and S. Wieczorek 



Is the astronomical forcing a reliable and unique pacemaker for Climate? 



28 



[Doedel et al, 2009] Doedel, E., Champneys, A., Der- 
cole, K, Fairgrieve, T., Kuznetsov, Y., Oldeman, B., 
Paffenroth, R., Sandstede, B., Wang, X., and Zhang, 
C. (2009). Auto: Software for continuation and 
bifurcation problems in ordinary differential equa- 
tions. Technical report, Montreal. 

[Donges et al., 2009] Donges, J. E, Zou, Y., Mar- 
wan, N., and Kurths, J. (2009). The backbone of 
the climate network. EPL (Europhysics Letters), 
87(4):48007. 

[Drazin and King, 1992] Drazin, P. G. and King, G. P., 
editors (1992). Interpretation of Time Series from 
Nonlinear Systems, volume 58. 

[Eckmann and Ruelle, 1985] Eckmann, J. P. and Ru- 
elle, D. (1985). Ergodic-theory of chaos and strange 
attractors. Reviews of Modem Physics, 57(3):617- 
656. 

[Ginoux and Rossetto, 2006] Ginoux, J.-M. and Ros- 
setto, B. (2006). Differential geometry and mechan- 
ics: Applications to chaotic dynamical systems. /. J. 
Bifurcation and Chaos, 16(4): 8 87-9 10. 

[Glass and Sun, 1994] Glass, L. and Sun, J. (1994). 
Periodic forcing of a limit-cycle oscillator: Eixed 
points, Arnold tongues, and the global organization 
of bifurcations. Phys. Rev. E, 50:5077-5084. 

[Glendinning and Wiersig, 1999] Glendinning, P. and 
Wiersig, J. (1999). Eine structure of mode-locked 
regions of the quasi-periodically forced circle map. 
Physics Letters A, 257(l-2):65-69. 

[Grasman et al., 2005] Grasman, J., Verhulst, E, and 
Shih, S. (2005). The Lyapunov exponents of the Van 
der Pol oscillator. Mathematical Methods in the Ap- 
plied Sciences, 28: 1131-1 139. 

[Grassberger and Procaccia, 1983] Grassberger, P. and 
Procaccia, I. (1983). Measuring the strangeness of 
strange attractors. Physica D: Nonlinear Phenom- 
ena, 9(1-2): 189-208. 

[Guckenheimer and Haiduc, 2005] Guckenheimer, J. 
and Haiduc, R. (2005). Canards at folded node. 
Mosc. Math. J, 5:91-103. 

[Guckenheimer et al., 2000] Guckenheimer, J., Hoff- 
man, K., and Weckesser, W. (2000). Numerical com- 
putation of canards. International Journal of Bifur- 
cation and Chaos, 10(1 2): 2269-2687. 



[Guckenheimer and Holmes, 1983] Guckenheimer, J. 
and Holmes, P. (1983). Nonlinear oscillations, dy- 
namical systems, and bifurcations of vector fields. 
Springer- Verlag, New York. 

[Haken, 1983] Haken, H. (1983). At least one lya- 
punov exponent vanishes if the trajectory of an at- 
tractor does not contain a fixed point. Physics Letters 
A, 94(2):71-72. 

[Hays et al., 1976] Hays, J. D., Imbrie, J., and Shack- 
leton, N. J. (1976). Variations in the Earth's orbit : 
Pacemaker of ice ages. Science, 194:1121-1132. 

[Hilborn, 2000] Hilborn, R. (2000). Chaos and Non- 
linear Dynamics: an Introduction for Scientists and 
Engineers. Oxford University Press. 

[Huybers, 2007] Huybers, R (2007). Glacial variability 
over the last two millions years: an extended depth- 
derived age model, continous obliquity pacing, and 
the Pleistocene progression. Quaternary Sci. Rev., 
26:37-55. 

[Imbrie and Imbrie, 1980] Imbrie, J. and Imbrie, J. Z. 
(1980). Modelling the climatic response to orbital 
variations. Science, 207:943-953. 

[Kantz and Schreiber, 2004] Kantz, H. and Schreiber, 
T. (2004). Nonlinear Time Series Analysis. Cam- 
bridge University Press, Cambridge, U.K., 2nd edi- 
tion. 

[Kaplan and Yorke, 1979] Kaplan, J. L. and Yorke, 
J. A. (1979). Chaotic behavior of multidimen- 
sional difference equations. In Functional differ- 
ential equations and approximation of fixed points 
(Proc. Summer School and Conf, Univ. Bonn, Bonn, 
1978), volume 730 of Lecture Notes in Math., pages 
204-227. Springer, Berlin. 

[Kloeden, 2000] Kloeden, R E. (2000). A Lyapunov 
function for puUback attractors of nonautonomous 
differential equations. Electronic J. Diff. Eqns, Conf. 
05, pages 91-102. 

[Kosmidis and Pakdaman, 2003] Kosmidis, E. K. and 
Pakdaman, K. (2003). An analysis of the reliability 
phenomenon in the fitzhugh-nagumo model. Journal 
of Computational Neuroscience, 14(l):5-22. 

[Langa et al., 2002] Langa, J. A., Robinson, J. C, and 
Suarez, A. (2002). Stability, instability, and bifur- 
cation phenomena in non- autonomous differential 
equations. Nonlinearity, 15(3): 1-17. 



Article submitted to Climate Dynamics 



B. De Saedeleer, M. Crucifix and S. Wieczorek 



Is the astronomical forcing a reliable and unique pacemaker for Climate? 



29 



[Laskar et al., 2004] Laskar, J., Robutel, P., Joutel, E, 
Boudin, R, Gastineau, M., Correia, A. C. M., and 
Levrard, B. (2004). A long-term numerical solution 
for the insolation quantities of the Earth. Astrom. 
Astropk, 42S:261-2S5. 

[Le Treut and Ghil, 1983] Le Trent, H. and Ghil, M. 
(1983). Orbital forcing, climatic interactions and 
glaciation cycles. J. Geophys. Res., 88(C9):5167- 
5190. 

[Lichtenberg and Lieberman, 1983] Lichtenberg, A. J. 
and Lieberman, M. A. (1983). Regular and stochas- 
tic motion. Springer- Verlag, New York. 

[Lisiecki and Raymo, 2005] Lisiecki, L. E. and 
Raymo, M. E. (2005). A pliocene-pleistocene stack 
of 57 globally distributed benthic 5^^0 records. 
Paleoceanography, 20( 1 ) . 

[Lisiecki and Raymo, 2007] Lisiecki, L. E. and 
Raymo, M. E. (2007). Plio-pleistocene climate 
evolution: trends and transitions in glacial cycle 
dynamics. Quaternary Science Reviews, 26(1- 
2):56-69. 

[Liu et al., 2005] Liu, H.-E, Dai, Z.-H., Li, W.-E, 
Gong, X., and Yu, Z.-H. (2005). Noise robust es- 
timates of the largest lyapunov exponent. Physics 
Letters A, 341(1-4): 119-127. 

[Luethi et al., 2008] Luethi, D., Le Eloch, M., Bere- 
iter, B., Blunier, T., Barnola, J.-M., Siegenthaler, U., 
Raynaud, D., Jouzel, J., Fischer, H., Kawamura, K., 
and Stocker, T. E. (2008). High-resolution carbon 
dioxide concentration record 650,000-800,000 years 
before present. Nature, 453(7193):379-382. 

[Marwan et al., 2009] Marwan, N., Donges, J. E, Zou, 
Y, Donner, R. V., and Kurths, J. (2009). Complex 
network approach for recurrence analysis of time se- 
ries. Physics Letters A, 373(46):4246-4254. 

[McCaffrey et al., 1992] McCaffrey, D. E, Ellner, S., 
Gallant, A. R., and Nychka, D. W. (Sep., 1992). Es- 
timating the lyapunov exponent of a chaotic system 
with nonparametric regression. Journal of the Amer- 
ican Statistical Association, 87(419):682-695. 

[Mettin et al., 1993] Mettin, R., Parlitz, U., and Lauter- 
born, W. (1993). Bifurcation structure of the driven 
van der pol oscillator. Int. J. Bifurcation & Chaos, 
(3): 1529-1555. 

Article submitted to Climate Dynamics 



[Milankovitch, 1941] Milankovitch, M. (1941). Kanon 
der Erdbestrahlung und Seine Anwendung auf das 
Eiszeitenproblem ( Canon of insolation and the ice- 
age problem). Konlishe Serbische Akademie, Bel- 
grad. 

[Miiller, 1995] Miiller, R C. (1995). Calculation of lya- 
punov exponents for dynamic systems with discon- 
tinuities. Chaos, Solitons & Fractals, 5(9): 1671- 
1681. 

[Oseledec, 1968] Oseledec, V. (1968). A multiplicative 
ergodic theorem: Ljapunov characteristic numbers 
for dynamical systems. Transactions of the Moscow 
Mathematical Society, 19:1 97-23 1 . 

[Osinga et al., 2000] Osinga, H., Wiersig, J., Glendin- 
ning. P., and Eeudel, U. (2000). Multi stability 
and nonsmooth bifurcations in the quasiperiodically 
forced circle map. ArXiv Nonlinear Sciences e- 
prints. 

[Ott, 2002] Ott, E. (2002). Chaos in Dynamical Sys- 
tems. Cambridge University Press. 

[Paillard, 1998] Paillard, D. (1998). The timing of 
pleistocene glaciations from a simple multiple- state 
climate model. Nature, 391:378-381. 

[Paillard and Parrenin, 2004] Paillard, D. and Parrenin, 
E. (2004). The Antarctic ice sheet and the triggering 
of deglaciations. Earth Planet. Sc. Lett, 227:263- 
271. 

[Parlitz and Lauterbom, 1987] Parlitz, U. and Lauter- 
born, W. (1987). Period-doubling cascades and 
devil's staircases of the driven van der pol oscillator. 
Physical Review A, 36(3). 

[Pikovsky et al., 2001] Pikovsky, A., Rosenblum, M., 
and Kurths, J. (2001). Synchronization A Universal 
Concept in Nonlinear Sciences. Cambridge Univer- 
sity Press, New York. 

[Rahmstorf etal., 2005] Rahmstorf, S., Crucifix, M., 
Ganopolski, A., Goosse, H., Kamenkovich, L, 
Knutti, R., Lohmann, G., Marsh, R., Mysak, L. A., 
Wang, Z., and Weaver, A. J. (2005). Thermoha- 
line circulation hysteresis: A model intercompari- 
son. Geophys. Res. Lett., 32(23). 

[Ramasubramanian and Sriram, 2000] 

Ramasubramanian, K. and Sriram, M. S. (2000). 

B. De Saedeleer, M. Crucifix and S. Wieczorek 



Is the astronomical forcing a reliable and unique pacemaker for Climate? 



30 



A comparative study of computation of lyapunov 
spectra with different algorithms. Physica D: 
Nonlinear Phenomena, 139(l-2):72-86. 

[Rosenstein et al., 1993] Rosenstein, M. T., Collins, 
J. J., and Luca, C. J. D. (1993). A practical 
method for calculating largest lyapunov exponents 
from small datasets. Physica D, (124). 

[Rossler, 1979] Rossler, O. E. (1979). An equation for 
hyperchaos. Physics Letters A, 71(2-3): 155-157. 

[Ruelle, 1990] Ruelle, D. (1990). Deterministic chaos: 
the science and the fiction. Proceedings of the Royal 
Society A, London, (427) :24 1-248. 

[Ruihong et al., 2008] Ruihong, L., Wei, X., and 
Shuang, L. (2008). Chaos control and synchroniza- 
tion of the 0^-van der pol system driven by exter- 
nal and parametric excitations. Nonlinear Dynamics, 
53(3):261-271. 

[Rulkov et al., 1995] Rulkov, N. E, Sushchik, M. M., 
Tsimring, L. S., and Abarbanel, H. D. I. (1995). Gen- 
eralized synchronization of chaos in directionally 
coupled chaotic systems. Phys. Rev. E, 51(2):980- 
994. 

[Saltzman, 2002] Saltzman, B. (2002). Dynamical Pa- 
leoclimatology: Generalized Theory of Global Cli- 
mate Change (International Geophysics). Academic 
Press. 

[Saltzman et al., 1984] Saltzman, B., Hansen, A. R., 
and Maasch, K. A. (1984). The late Quater- 
nary glaciations as the response of a 3 -component 
feedback- system to Earth-orbital forcing. Journal of 
the Atmospheric Sciences, 41(23):3380-3389. 

[Saltzman and Maasch, 1990] Saltzman, B. and 
Maasch, K. A. (1990). A first-order global model 
of late Cenozoic climate. Trans. R. Soc. Edinburgh 
Earth Sci,Sl:315-325. 

[Saltzman and Maasch, 1991] Saltzman, B. and 
Maasch, K. A. (1991). A first-order global model of 
late Cenozoic climate. II further analysis based on 
a simplification of the CO2 dynamics. Clim. Dyn., 
5:201-210. 

[Savi, 2005] Savi, M. A. (2005). Chaos and order 
in biomedical rhythms. Journal of the Brazilian 
Society of Mechanical Sciences and Engineering, 
27(2): 157-169. 

Article submitted to Climate Dynamics 



[Shimada and Nagashima, 1979] Shimada, I. and Na- 
gashima, T. (1979). A numerical approach to ergodic 
problem of dissipative dynamical systems. Progr 
Theoret. Phys., 61(6): 1605-1616. 

[Strogatz, 1994] Strogatz, S. H. (1994). Nonlinear Dy- 
namics And Chaos: With Applications To Physics, 
Biology, Chemistry, And Engineering (Studies in 
Nonlinearity). Studies in nonlinearity. Perseus 
Books Group, 1 edition. 

[Svensson and Coombes, 2009] Svensson, C.-M. and 
Coombes, S. (2009). Mode locking in a spatially 
extended neuron model: active soma and compart- 
mental tree. International Journal of Bifurcation and 
Chaos, 19(8):2597-2607. 

[Theiler, 1990] Theiler, J. (1990). Estimating fractal 
dimension. /. Opt. Soc. Am. A, 7(6): 1055-1073. 

[Tsiganis et al., 1999] Tsiganis, K., Anastasiadis, A., 
and Varvoglis, H. (1999). Effective lyapunov num- 
bers and correlation dimensions in a 3-D hamiltonian 
system. In Henrard, J. and Eerraz-Mello, S., editors, 
lAU -Colloquium No 172 - The impact of modern dy- 
namics in Astronomy, pages 447-448. 

[Tziperman and Gildor, 2003] Tziperman, E. and 
Gildor, H. (2003). On the mid-Pleistocene transtion 
to 100-kyr glacial cycles and the asymmetry between 
glaciation and deglaciation times. Paleoceanogra- 
phy, 18(1):1001. 

[Tziperman et al., 2006] Tziperman, E., Raymo, M. E., 
Huybers, P., and Wunsch, C. (2006). Consequences 
of pacing the Pleistocene 100 kyr ice ages by non- 
linear phase locking to Milankovitch forcing. Paleo- 
ceanography, 21:PA4206. 

[van der Pol, 1926] van der Pol, B. (1926). On relax- 
ation oscillations. Phil. Mag., 2(ll):978-992. 

[Wieczorek, 2009] Wieczorek, S. (2009). Stochastic 
bifurcation in noise-driven lasers and hopf oscilla- 
tors. Phys. Rev. E, 79(3):036209. 

[Wieczorek, 2011] Wieczorek, S. M. (2011). Noise 
synchronisation and stochastic bifurcations in lasers. 
http.V/arxiv org/abs/1 1 04. 4052 . 

[Wiggins, 2003] Wiggins, S. (2003). Introduction to 
Applied Nonlinear Dynamical Systems and Chaos. 
Texts in Applied Mathematics. Springer, 2nd edition. 

B. De Saedeleer, M. Crucifix and S. Wieczorek 



Is the astronomical forcing a reliable and unique pacemaker for Climate? 



31 



[Wolf etal., 1985] Wolf, A., Swift, J. B., Swinney, 
H. L., and Vastano, J. A. (1985). Determining lya- 
punov exponents from a time series. Physica D: 
Nonlinear Phenomena, 16(3):285-317. 

[Wu et al., 2006] Wu, L., Zhu, S., and Li, J. (2006). 
Synchronization on fast and slow dynamics in drive- 
response systems. Physica D: Nonlinear Phenom- 
ena, 223(2):208-213. 



Article submitted to Climate Dynamics 



B. De Saedeleer, M. Crucifix and S. Wieczorek 



